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