001:       SUBROUTINE ZLALSD( 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:       DOUBLE PRECISION   RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
016:       COMPLEX*16         B( LDB, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  ZLALSD 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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*16 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) DOUBLE PRECISION
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*16 array, dimension at least
088: *         (N * NRHS).
089: *
090: *  RWORK  (workspace) DOUBLE PRECISION 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 at least
096: *         (3*N*NLVL + 11*N).
097: *
098: *  INFO   (output) INTEGER
099: *         = 0:  successful exit.
100: *         < 0:  if INFO = -i, the i-th argument had an illegal value.
101: *         > 0:  The algorithm failed to compute an singular value while
102: *               working on the submatrix lying in rows and columns
103: *               INFO/(N+1) through MOD(INFO,N+1).
104: *
105: *  Further Details
106: *  ===============
107: *
108: *  Based on contributions by
109: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
110: *       California at Berkeley, USA
111: *     Osni Marques, LBNL/NERSC, USA
112: *
113: *  =====================================================================
114: *
115: *     .. Parameters ..
116:       DOUBLE PRECISION   ZERO, ONE, TWO
117:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
118:       COMPLEX*16         CZERO
119:       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
120: *     ..
121: *     .. Local Scalars ..
122:       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
123:      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
124:      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
125:      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
126:      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
127:      $                   U, VT, Z
128:       DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
129: *     ..
130: *     .. External Functions ..
131:       INTEGER            IDAMAX
132:       DOUBLE PRECISION   DLAMCH, DLANST
133:       EXTERNAL           IDAMAX, DLAMCH, DLANST
134: *     ..
135: *     .. External Subroutines ..
136:       EXTERNAL           DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
137:      $                   DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
138:      $                   ZLASCL, ZLASET
139: *     ..
140: *     .. Intrinsic Functions ..
141:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
142: *     ..
143: *     .. Executable Statements ..
144: *
145: *     Test the input parameters.
146: *
147:       INFO = 0
148: *
149:       IF( N.LT.0 ) THEN
150:          INFO = -3
151:       ELSE IF( NRHS.LT.1 ) THEN
152:          INFO = -4
153:       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
154:          INFO = -8
155:       END IF
156:       IF( INFO.NE.0 ) THEN
157:          CALL XERBLA( 'ZLALSD', -INFO )
158:          RETURN
159:       END IF
160: *
161:       EPS = DLAMCH( 'Epsilon' )
162: *
163: *     Set up the tolerance.
164: *
165:       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
166:          RCND = EPS
167:       ELSE
168:          RCND = RCOND
169:       END IF
170: *
171:       RANK = 0
172: *
173: *     Quick return if possible.
174: *
175:       IF( N.EQ.0 ) THEN
176:          RETURN
177:       ELSE IF( N.EQ.1 ) THEN
178:          IF( D( 1 ).EQ.ZERO ) THEN
179:             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
180:          ELSE
181:             RANK = 1
182:             CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
183:             D( 1 ) = ABS( D( 1 ) )
184:          END IF
185:          RETURN
186:       END IF
187: *
188: *     Rotate the matrix if it is lower bidiagonal.
189: *
190:       IF( UPLO.EQ.'L' ) THEN
191:          DO 10 I = 1, N - 1
192:             CALL DLARTG( D( I ), E( I ), CS, SN, R )
193:             D( I ) = R
194:             E( I ) = SN*D( I+1 )
195:             D( I+1 ) = CS*D( I+1 )
196:             IF( NRHS.EQ.1 ) THEN
197:                CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
198:             ELSE
199:                RWORK( I*2-1 ) = CS
200:                RWORK( I*2 ) = SN
201:             END IF
202:    10    CONTINUE
203:          IF( NRHS.GT.1 ) THEN
204:             DO 30 I = 1, NRHS
205:                DO 20 J = 1, N - 1
206:                   CS = RWORK( J*2-1 )
207:                   SN = RWORK( J*2 )
208:                   CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
209:    20          CONTINUE
210:    30       CONTINUE
211:          END IF
212:       END IF
213: *
214: *     Scale.
215: *
216:       NM1 = N - 1
217:       ORGNRM = DLANST( 'M', N, D, E )
218:       IF( ORGNRM.EQ.ZERO ) THEN
219:          CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
220:          RETURN
221:       END IF
222: *
223:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
224:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
225: *
226: *     If N is smaller than the minimum divide size SMLSIZ, then solve
227: *     the problem with another solver.
228: *
229:       IF( N.LE.SMLSIZ ) THEN
230:          IRWU = 1
231:          IRWVT = IRWU + N*N
232:          IRWWRK = IRWVT + N*N
233:          IRWRB = IRWWRK
234:          IRWIB = IRWRB + N*NRHS
235:          IRWB = IRWIB + N*NRHS
236:          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
237:          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
238:          CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
239:      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
240:      $                RWORK( IRWWRK ), INFO )
241:          IF( INFO.NE.0 ) THEN
242:             RETURN
243:          END IF
244: *
245: *        In the real version, B is passed to DLASDQ and multiplied
246: *        internally by Q'. Here B is complex and that product is
247: *        computed below in two steps (real and imaginary parts).
248: *
249:          J = IRWB - 1
250:          DO 50 JCOL = 1, NRHS
251:             DO 40 JROW = 1, N
252:                J = J + 1
253:                RWORK( J ) = DBLE( B( JROW, JCOL ) )
254:    40       CONTINUE
255:    50    CONTINUE
256:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
257:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
258:          J = IRWB - 1
259:          DO 70 JCOL = 1, NRHS
260:             DO 60 JROW = 1, N
261:                J = J + 1
262:                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
263:    60       CONTINUE
264:    70    CONTINUE
265:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
266:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
267:          JREAL = IRWRB - 1
268:          JIMAG = IRWIB - 1
269:          DO 90 JCOL = 1, NRHS
270:             DO 80 JROW = 1, N
271:                JREAL = JREAL + 1
272:                JIMAG = JIMAG + 1
273:                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
274:      $                           RWORK( JIMAG ) )
275:    80       CONTINUE
276:    90    CONTINUE
277: *
278:          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
279:          DO 100 I = 1, N
280:             IF( D( I ).LE.TOL ) THEN
281:                CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
282:             ELSE
283:                CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
284:      $                      LDB, INFO )
285:                RANK = RANK + 1
286:             END IF
287:   100    CONTINUE
288: *
289: *        Since B is complex, the following call to DGEMM is performed
290: *        in two steps (real and imaginary parts). That is for V * B
291: *        (in the real version of the code V' is stored in WORK).
292: *
293: *        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
294: *    $               WORK( NWORK ), N )
295: *
296:          J = IRWB - 1
297:          DO 120 JCOL = 1, NRHS
298:             DO 110 JROW = 1, N
299:                J = J + 1
300:                RWORK( J ) = DBLE( B( JROW, JCOL ) )
301:   110       CONTINUE
302:   120    CONTINUE
303:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
304:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
305:          J = IRWB - 1
306:          DO 140 JCOL = 1, NRHS
307:             DO 130 JROW = 1, N
308:                J = J + 1
309:                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
310:   130       CONTINUE
311:   140    CONTINUE
312:          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
313:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
314:          JREAL = IRWRB - 1
315:          JIMAG = IRWIB - 1
316:          DO 160 JCOL = 1, NRHS
317:             DO 150 JROW = 1, N
318:                JREAL = JREAL + 1
319:                JIMAG = JIMAG + 1
320:                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
321:      $                           RWORK( JIMAG ) )
322:   150       CONTINUE
323:   160    CONTINUE
324: *
325: *        Unscale.
326: *
327:          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
328:          CALL DLASRT( 'D', N, D, INFO )
329:          CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
330: *
331:          RETURN
332:       END IF
333: *
334: *     Book-keeping and setting up some constants.
335: *
336:       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
337: *
338:       SMLSZP = SMLSIZ + 1
339: *
340:       U = 1
341:       VT = 1 + SMLSIZ*N
342:       DIFL = VT + SMLSZP*N
343:       DIFR = DIFL + NLVL*N
344:       Z = DIFR + NLVL*N*2
345:       C = Z + NLVL*N
346:       S = C + N
347:       POLES = S + N
348:       GIVNUM = POLES + 2*NLVL*N
349:       NRWORK = GIVNUM + 2*NLVL*N
350:       BX = 1
351: *
352:       IRWRB = NRWORK
353:       IRWIB = IRWRB + SMLSIZ*NRHS
354:       IRWB = IRWIB + SMLSIZ*NRHS
355: *
356:       SIZEI = 1 + N
357:       K = SIZEI + N
358:       GIVPTR = K + N
359:       PERM = GIVPTR + N
360:       GIVCOL = PERM + NLVL*N
361:       IWK = GIVCOL + NLVL*N*2
362: *
363:       ST = 1
364:       SQRE = 0
365:       ICMPQ1 = 1
366:       ICMPQ2 = 0
367:       NSUB = 0
368: *
369:       DO 170 I = 1, N
370:          IF( ABS( D( I ) ).LT.EPS ) THEN
371:             D( I ) = SIGN( EPS, D( I ) )
372:          END IF
373:   170 CONTINUE
374: *
375:       DO 240 I = 1, NM1
376:          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
377:             NSUB = NSUB + 1
378:             IWORK( NSUB ) = ST
379: *
380: *           Subproblem found. First determine its size and then
381: *           apply divide and conquer on it.
382: *
383:             IF( I.LT.NM1 ) THEN
384: *
385: *              A subproblem with E(I) small for I < NM1.
386: *
387:                NSIZE = I - ST + 1
388:                IWORK( SIZEI+NSUB-1 ) = NSIZE
389:             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
390: *
391: *              A subproblem with E(NM1) not too small but I = NM1.
392: *
393:                NSIZE = N - ST + 1
394:                IWORK( SIZEI+NSUB-1 ) = NSIZE
395:             ELSE
396: *
397: *              A subproblem with E(NM1) small. This implies an
398: *              1-by-1 subproblem at D(N), which is not solved
399: *              explicitly.
400: *
401:                NSIZE = I - ST + 1
402:                IWORK( SIZEI+NSUB-1 ) = NSIZE
403:                NSUB = NSUB + 1
404:                IWORK( NSUB ) = N
405:                IWORK( SIZEI+NSUB-1 ) = 1
406:                CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
407:             END IF
408:             ST1 = ST - 1
409:             IF( NSIZE.EQ.1 ) THEN
410: *
411: *              This is a 1-by-1 subproblem and is not solved
412: *              explicitly.
413: *
414:                CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
415:             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
416: *
417: *              This is a small subproblem and is solved by DLASDQ.
418: *
419:                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
420:      $                      RWORK( VT+ST1 ), N )
421:                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
422:      $                      RWORK( U+ST1 ), N )
423:                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
424:      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
425:      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
426:      $                      INFO )
427:                IF( INFO.NE.0 ) THEN
428:                   RETURN
429:                END IF
430: *
431: *              In the real version, B is passed to DLASDQ and multiplied
432: *              internally by Q'. Here B is complex and that product is
433: *              computed below in two steps (real and imaginary parts).
434: *
435:                J = IRWB - 1
436:                DO 190 JCOL = 1, NRHS
437:                   DO 180 JROW = ST, ST + NSIZE - 1
438:                      J = J + 1
439:                      RWORK( J ) = DBLE( B( JROW, JCOL ) )
440:   180             CONTINUE
441:   190          CONTINUE
442:                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
443:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
444:      $                     ZERO, RWORK( IRWRB ), NSIZE )
445:                J = IRWB - 1
446:                DO 210 JCOL = 1, NRHS
447:                   DO 200 JROW = ST, ST + NSIZE - 1
448:                      J = J + 1
449:                      RWORK( J ) = DIMAG( B( JROW, JCOL ) )
450:   200             CONTINUE
451:   210          CONTINUE
452:                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
453:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
454:      $                     ZERO, RWORK( IRWIB ), NSIZE )
455:                JREAL = IRWRB - 1
456:                JIMAG = IRWIB - 1
457:                DO 230 JCOL = 1, NRHS
458:                   DO 220 JROW = ST, ST + NSIZE - 1
459:                      JREAL = JREAL + 1
460:                      JIMAG = JIMAG + 1
461:                      B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
462:      $                                 RWORK( JIMAG ) )
463:   220             CONTINUE
464:   230          CONTINUE
465: *
466:                CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
467:      $                      WORK( BX+ST1 ), N )
468:             ELSE
469: *
470: *              A large problem. Solve it using divide and conquer.
471: *
472:                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
473:      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
474:      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
475:      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
476:      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
477:      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
478:      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
479:      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
480:      $                      IWORK( IWK ), INFO )
481:                IF( INFO.NE.0 ) THEN
482:                   RETURN
483:                END IF
484:                BXST = BX + ST1
485:                CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
486:      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
487:      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
488:      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
489:      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
490:      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
491:      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
492:      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
493:      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
494:                IF( INFO.NE.0 ) THEN
495:                   RETURN
496:                END IF
497:             END IF
498:             ST = I + 1
499:          END IF
500:   240 CONTINUE
501: *
502: *     Apply the singular values and treat the tiny ones as zero.
503: *
504:       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
505: *
506:       DO 250 I = 1, N
507: *
508: *        Some of the elements in D can be negative because 1-by-1
509: *        subproblems were not solved explicitly.
510: *
511:          IF( ABS( D( I ) ).LE.TOL ) THEN
512:             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
513:          ELSE
514:             RANK = RANK + 1
515:             CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
516:      $                   WORK( BX+I-1 ), N, INFO )
517:          END IF
518:          D( I ) = ABS( D( I ) )
519:   250 CONTINUE
520: *
521: *     Now apply back the right singular vectors.
522: *
523:       ICMPQ2 = 1
524:       DO 320 I = 1, NSUB
525:          ST = IWORK( I )
526:          ST1 = ST - 1
527:          NSIZE = IWORK( SIZEI+I-1 )
528:          BXST = BX + ST1
529:          IF( NSIZE.EQ.1 ) THEN
530:             CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
531:          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
532: *
533: *           Since B and BX are complex, the following call to DGEMM
534: *           is performed in two steps (real and imaginary parts).
535: *
536: *           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
537: *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
538: *    $                  B( ST, 1 ), LDB )
539: *
540:             J = BXST - N - 1
541:             JREAL = IRWB - 1
542:             DO 270 JCOL = 1, NRHS
543:                J = J + N
544:                DO 260 JROW = 1, NSIZE
545:                   JREAL = JREAL + 1
546:                   RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
547:   260          CONTINUE
548:   270       CONTINUE
549:             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
550:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
551:      $                  RWORK( IRWRB ), NSIZE )
552:             J = BXST - N - 1
553:             JIMAG = IRWB - 1
554:             DO 290 JCOL = 1, NRHS
555:                J = J + N
556:                DO 280 JROW = 1, NSIZE
557:                   JIMAG = JIMAG + 1
558:                   RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
559:   280          CONTINUE
560:   290       CONTINUE
561:             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
562:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
563:      $                  RWORK( IRWIB ), NSIZE )
564:             JREAL = IRWRB - 1
565:             JIMAG = IRWIB - 1
566:             DO 310 JCOL = 1, NRHS
567:                DO 300 JROW = ST, ST + NSIZE - 1
568:                   JREAL = JREAL + 1
569:                   JIMAG = JIMAG + 1
570:                   B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
571:      $                              RWORK( JIMAG ) )
572:   300          CONTINUE
573:   310       CONTINUE
574:          ELSE
575:             CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
576:      $                   B( ST, 1 ), LDB, 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:   320 CONTINUE
589: *
590: *     Unscale and sort the singular values.
591: *
592:       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
593:       CALL DLASRT( 'D', N, D, INFO )
594:       CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
595: *
596:       RETURN
597: *
598: *     End of ZLALSD
599: *
600:       END
601: