001:       SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
002:      $                   WORK, LWORK, RWORK, IWORK, INFO )
003: *
004: *  -- LAPACK driver 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:       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
011:       DOUBLE PRECISION   RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       DOUBLE PRECISION   RWORK( * ), S( * )
016:       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  ZGELSD computes the minimum-norm solution to a real linear least
023: *  squares problem:
024: *      minimize 2-norm(| b - A*x |)
025: *  using the singular value decomposition (SVD) of A. A is an M-by-N
026: *  matrix which may be rank-deficient.
027: *
028: *  Several right hand side vectors b and solution vectors x can be
029: *  handled in a single call; they are stored as the columns of the
030: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
031: *  matrix X.
032: *
033: *  The problem is solved in three steps:
034: *  (1) Reduce the coefficient matrix A to bidiagonal form with
035: *      Householder tranformations, reducing the original problem
036: *      into a "bidiagonal least squares problem" (BLS)
037: *  (2) Solve the BLS using a divide and conquer approach.
038: *  (3) Apply back all the Householder tranformations to solve
039: *      the original least squares problem.
040: *
041: *  The effective rank of A is determined by treating as zero those
042: *  singular values which are less than RCOND times the largest singular
043: *  value.
044: *
045: *  The divide and conquer algorithm makes very mild assumptions about
046: *  floating point arithmetic. It will work on machines with a guard
047: *  digit in add/subtract, or on those binary machines without guard
048: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
049: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
050: *  without guard digits, but we know of none.
051: *
052: *  Arguments
053: *  =========
054: *
055: *  M       (input) INTEGER
056: *          The number of rows of the matrix A. M >= 0.
057: *
058: *  N       (input) INTEGER
059: *          The number of columns of the matrix A. N >= 0.
060: *
061: *  NRHS    (input) INTEGER
062: *          The number of right hand sides, i.e., the number of columns
063: *          of the matrices B and X. NRHS >= 0.
064: *
065: *  A       (input) COMPLEX*16 array, dimension (LDA,N)
066: *          On entry, the M-by-N matrix A.
067: *          On exit, A has been destroyed.
068: *
069: *  LDA     (input) INTEGER
070: *          The leading dimension of the array A. LDA >= max(1,M).
071: *
072: *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
073: *          On entry, the M-by-NRHS right hand side matrix B.
074: *          On exit, B is overwritten by the N-by-NRHS solution matrix X.
075: *          If m >= n and RANK = n, the residual sum-of-squares for
076: *          the solution in the i-th column is given by the sum of
077: *          squares of the modulus of elements n+1:m in that column.
078: *
079: *  LDB     (input) INTEGER
080: *          The leading dimension of the array B.  LDB >= max(1,M,N).
081: *
082: *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
083: *          The singular values of A in decreasing order.
084: *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
085: *
086: *  RCOND   (input) DOUBLE PRECISION
087: *          RCOND is used to determine the effective rank of A.
088: *          Singular values S(i) <= RCOND*S(1) are treated as zero.
089: *          If RCOND < 0, machine precision is used instead.
090: *
091: *  RANK    (output) INTEGER
092: *          The effective rank of A, i.e., the number of singular values
093: *          which are greater than RCOND*S(1).
094: *
095: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
096: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
097: *
098: *  LWORK   (input) INTEGER
099: *          The dimension of the array WORK. LWORK must be at least 1.
100: *          The exact minimum amount of workspace needed depends on M,
101: *          N and NRHS. As long as LWORK is at least
102: *              2*N + N*NRHS
103: *          if M is greater than or equal to N or
104: *              2*M + M*NRHS
105: *          if M is less than N, the code will execute correctly.
106: *          For good performance, LWORK should generally be larger.
107: *
108: *          If LWORK = -1, then a workspace query is assumed; the routine
109: *          only calculates the optimal size of the array WORK and the
110: *          minimum sizes of the arrays RWORK and IWORK, and returns
111: *          these values as the first entries of the WORK, RWORK and
112: *          IWORK arrays, and no error message related to LWORK is issued
113: *          by XERBLA.
114: *
115: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
116: *          LRWORK >=
117: *              10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
118: *             (SMLSIZ+1)**2
119: *          if M is greater than or equal to N or
120: *             10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
121: *             (SMLSIZ+1)**2
122: *          if M is less than N, the code will execute correctly.
123: *          SMLSIZ is returned by ILAENV and is equal to the maximum
124: *          size of the subproblems at the bottom of the computation
125: *          tree (usually about 25), and
126: *             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
127: *          On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
128: *
129: *  IWORK   (workspace) INTEGER array, dimension (MAX(1,LIWORK))
130: *          LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
131: *          where MINMN = MIN( M,N ).
132: *          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
133: *
134: *  INFO    (output) INTEGER
135: *          = 0: successful exit
136: *          < 0: if INFO = -i, the i-th argument had an illegal value.
137: *          > 0:  the algorithm for computing the SVD failed to converge;
138: *                if INFO = i, i off-diagonal elements of an intermediate
139: *                bidiagonal form did not converge to zero.
140: *
141: *  Further Details
142: *  ===============
143: *
144: *  Based on contributions by
145: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
146: *       California at Berkeley, USA
147: *     Osni Marques, LBNL/NERSC, USA
148: *
149: *  =====================================================================
150: *
151: *     .. Parameters ..
152:       DOUBLE PRECISION   ZERO, ONE, TWO
153:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
154:       COMPLEX*16         CZERO
155:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
156: *     ..
157: *     .. Local Scalars ..
158:       LOGICAL            LQUERY
159:       INTEGER            IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
160:      $                   LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
161:      $                   MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
162:       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
163: *     ..
164: *     .. External Subroutines ..
165:       EXTERNAL           DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF,
166:      $                   ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR,
167:      $                   ZUNMLQ, ZUNMQR
168: *     ..
169: *     .. External Functions ..
170:       INTEGER            ILAENV
171:       DOUBLE PRECISION   DLAMCH, ZLANGE
172:       EXTERNAL           ILAENV, DLAMCH, ZLANGE
173: *     ..
174: *     .. Intrinsic Functions ..
175:       INTRINSIC          INT, LOG, MAX, MIN, DBLE
176: *     ..
177: *     .. Executable Statements ..
178: *
179: *     Test the input arguments.
180: *
181:       INFO = 0
182:       MINMN = MIN( M, N )
183:       MAXMN = MAX( M, N )
184:       LQUERY = ( LWORK.EQ.-1 )
185:       IF( M.LT.0 ) THEN
186:          INFO = -1
187:       ELSE IF( N.LT.0 ) THEN
188:          INFO = -2
189:       ELSE IF( NRHS.LT.0 ) THEN
190:          INFO = -3
191:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
192:          INFO = -5
193:       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
194:          INFO = -7
195:       END IF
196: *
197: *     Compute workspace.
198: *     (Note: Comments in the code beginning "Workspace:" describe the
199: *     minimal amount of workspace needed at that point in the code,
200: *     as well as the preferred amount for good performance.
201: *     NB refers to the optimal block size for the immediately
202: *     following subroutine, as returned by ILAENV.)
203: *
204:       IF( INFO.EQ.0 ) THEN
205:          MINWRK = 1
206:          MAXWRK = 1
207:          LIWORK = 1
208:          LRWORK = 1
209:          IF( MINMN.GT.0 ) THEN
210:             SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 )
211:             MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 )
212:             NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) /
213:      $                  LOG( TWO ) ) + 1, 0 )
214:             LIWORK = 3*MINMN*NLVL + 11*MINMN
215:             MM = M
216:             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
217: *
218: *              Path 1a - overdetermined, with many more rows than
219: *                        columns.
220: *
221:                MM = N
222:                MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N,
223:      $                       -1, -1 ) )
224:                MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M,
225:      $                       NRHS, N, -1 ) )
226:             END IF
227:             IF( M.GE.N ) THEN
228: *
229: *              Path 1 - overdetermined or exactly determined.
230: *
231:                LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
232:      $                  ( SMLSIZ + 1 )**2
233:                MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
234:      $                       'ZGEBRD', ' ', MM, N, -1, -1 ) )
235:                MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
236:      $                       'QLC', MM, NRHS, N, -1 ) )
237:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
238:      $                       'ZUNMBR', 'PLN', N, NRHS, N, -1 ) )
239:                MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
240:                MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
241:             END IF
242:             IF( N.GT.M ) THEN
243:                LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
244:      $                  ( SMLSIZ + 1 )**2
245:                IF( N.GE.MNTHR ) THEN
246: *
247: *                 Path 2a - underdetermined, with many more columns
248: *                           than rows.
249: *
250:                   MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
251:      $                     -1 )
252:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
253:      $                          'ZGEBRD', ' ', M, M, -1, -1 ) )
254:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
255:      $                          'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
256:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
257:      $                          'ZUNMLQ', 'LC', N, NRHS, M, -1 ) )
258:                   IF( NRHS.GT.1 ) THEN
259:                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
260:                   ELSE
261:                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
262:                   END IF
263:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
264: !     XXX: Ensure the Path 2a case below is triggered.  The workspace
265: !     calculation should use queries for all routines eventually.
266:                   MAXWRK = MAX( MAXWRK,
267:      $                 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
268:                ELSE
269: *
270: *                 Path 2 - underdetermined.
271: *
272:                   MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
273:      $                     N, -1, -1 )
274:                   MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
275:      $                          'QLC', M, NRHS, M, -1 ) )
276:                   MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR',
277:      $                          'PLN', N, NRHS, M, -1 ) )
278:                   MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
279:                END IF
280:                MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
281:             END IF
282:          END IF
283:          MINWRK = MIN( MINWRK, MAXWRK )
284:          WORK( 1 ) = MAXWRK
285:          IWORK( 1 ) = LIWORK
286:          RWORK( 1 ) = LRWORK
287: *
288:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
289:             INFO = -12
290:          END IF
291:       END IF
292: *
293:       IF( INFO.NE.0 ) THEN
294:          CALL XERBLA( 'ZGELSD', -INFO )
295:          RETURN
296:       ELSE IF( LQUERY ) THEN
297:          RETURN
298:       END IF
299: *
300: *     Quick return if possible.
301: *
302:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
303:          RANK = 0
304:          RETURN
305:       END IF
306: *
307: *     Get machine parameters.
308: *
309:       EPS = DLAMCH( 'P' )
310:       SFMIN = DLAMCH( 'S' )
311:       SMLNUM = SFMIN / EPS
312:       BIGNUM = ONE / SMLNUM
313:       CALL DLABAD( SMLNUM, BIGNUM )
314: *
315: *     Scale A if max entry outside range [SMLNUM,BIGNUM].
316: *
317:       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
318:       IASCL = 0
319:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
320: *
321: *        Scale matrix norm up to SMLNUM
322: *
323:          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
324:          IASCL = 1
325:       ELSE IF( ANRM.GT.BIGNUM ) THEN
326: *
327: *        Scale matrix norm down to BIGNUM.
328: *
329:          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
330:          IASCL = 2
331:       ELSE IF( ANRM.EQ.ZERO ) THEN
332: *
333: *        Matrix all zero. Return zero solution.
334: *
335:          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
336:          CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
337:          RANK = 0
338:          GO TO 10
339:       END IF
340: *
341: *     Scale B if max entry outside range [SMLNUM,BIGNUM].
342: *
343:       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
344:       IBSCL = 0
345:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
346: *
347: *        Scale matrix norm up to SMLNUM.
348: *
349:          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
350:          IBSCL = 1
351:       ELSE IF( BNRM.GT.BIGNUM ) THEN
352: *
353: *        Scale matrix norm down to BIGNUM.
354: *
355:          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
356:          IBSCL = 2
357:       END IF
358: *
359: *     If M < N make sure B(M+1:N,:) = 0
360: *
361:       IF( M.LT.N )
362:      $   CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
363: *
364: *     Overdetermined case.
365: *
366:       IF( M.GE.N ) THEN
367: *
368: *        Path 1 - overdetermined or exactly determined.
369: *
370:          MM = M
371:          IF( M.GE.MNTHR ) THEN
372: *
373: *           Path 1a - overdetermined, with many more rows than columns
374: *
375:             MM = N
376:             ITAU = 1
377:             NWORK = ITAU + N
378: *
379: *           Compute A=Q*R.
380: *           (RWorkspace: need N)
381: *           (CWorkspace: need N, prefer N*NB)
382: *
383:             CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
384:      $                   LWORK-NWORK+1, INFO )
385: *
386: *           Multiply B by transpose(Q).
387: *           (RWorkspace: need N)
388: *           (CWorkspace: need NRHS, prefer NRHS*NB)
389: *
390:             CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
391:      $                   LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
392: *
393: *           Zero out below R.
394: *
395:             IF( N.GT.1 ) THEN
396:                CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
397:      $                      LDA )
398:             END IF
399:          END IF
400: *
401:          ITAUQ = 1
402:          ITAUP = ITAUQ + N
403:          NWORK = ITAUP + N
404:          IE = 1
405:          NRWORK = IE + N
406: *
407: *        Bidiagonalize R in A.
408: *        (RWorkspace: need N)
409: *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
410: *
411:          CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
412:      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
413:      $                INFO )
414: *
415: *        Multiply B by transpose of left bidiagonalizing vectors of R.
416: *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
417: *
418:          CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
419:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
420: *
421: *        Solve the bidiagonal least squares problem.
422: *
423:          CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
424:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
425:      $                IWORK, INFO )
426:          IF( INFO.NE.0 ) THEN
427:             GO TO 10
428:          END IF
429: *
430: *        Multiply B by right bidiagonalizing vectors of R.
431: *
432:          CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
433:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
434: *
435:       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
436:      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
437: *
438: *        Path 2a - underdetermined, with many more columns than rows
439: *        and sufficient workspace for an efficient algorithm.
440: *
441:          LDWORK = M
442:          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
443:      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
444:          ITAU = 1
445:          NWORK = M + 1
446: *
447: *        Compute A=L*Q.
448: *        (CWorkspace: need 2*M, prefer M+M*NB)
449: *
450:          CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
451:      $                LWORK-NWORK+1, INFO )
452:          IL = NWORK
453: *
454: *        Copy L to WORK(IL), zeroing out above its diagonal.
455: *
456:          CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
457:          CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
458:      $                LDWORK )
459:          ITAUQ = IL + LDWORK*M
460:          ITAUP = ITAUQ + M
461:          NWORK = ITAUP + M
462:          IE = 1
463:          NRWORK = IE + M
464: *
465: *        Bidiagonalize L in WORK(IL).
466: *        (RWorkspace: need M)
467: *        (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
468: *
469:          CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
470:      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
471:      $                LWORK-NWORK+1, INFO )
472: *
473: *        Multiply B by transpose of left bidiagonalizing vectors of L.
474: *        (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
475: *
476:          CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
477:      $                WORK( ITAUQ ), B, LDB, WORK( NWORK ),
478:      $                LWORK-NWORK+1, INFO )
479: *
480: *        Solve the bidiagonal least squares problem.
481: *
482:          CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
483:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
484:      $                IWORK, INFO )
485:          IF( INFO.NE.0 ) THEN
486:             GO TO 10
487:          END IF
488: *
489: *        Multiply B by right bidiagonalizing vectors of L.
490: *
491:          CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
492:      $                WORK( ITAUP ), B, LDB, WORK( NWORK ),
493:      $                LWORK-NWORK+1, INFO )
494: *
495: *        Zero out below first M rows of B.
496: *
497:          CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
498:          NWORK = ITAU + M
499: *
500: *        Multiply transpose(Q) by B.
501: *        (CWorkspace: need NRHS, prefer NRHS*NB)
502: *
503:          CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
504:      $                LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
505: *
506:       ELSE
507: *
508: *        Path 2 - remaining underdetermined cases.
509: *
510:          ITAUQ = 1
511:          ITAUP = ITAUQ + M
512:          NWORK = ITAUP + M
513:          IE = 1
514:          NRWORK = IE + M
515: *
516: *        Bidiagonalize A.
517: *        (RWorkspace: need M)
518: *        (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
519: *
520:          CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
521:      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
522:      $                INFO )
523: *
524: *        Multiply B by transpose of left bidiagonalizing vectors.
525: *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
526: *
527:          CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
528:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
529: *
530: *        Solve the bidiagonal least squares problem.
531: *
532:          CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
533:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
534:      $                IWORK, INFO )
535:          IF( INFO.NE.0 ) THEN
536:             GO TO 10
537:          END IF
538: *
539: *        Multiply B by right bidiagonalizing vectors of A.
540: *
541:          CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
542:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
543: *
544:       END IF
545: *
546: *     Undo scaling.
547: *
548:       IF( IASCL.EQ.1 ) THEN
549:          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
550:          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
551:      $                INFO )
552:       ELSE IF( IASCL.EQ.2 ) THEN
553:          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
554:          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
555:      $                INFO )
556:       END IF
557:       IF( IBSCL.EQ.1 ) THEN
558:          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
559:       ELSE IF( IBSCL.EQ.2 ) THEN
560:          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
561:       END IF
562: *
563:    10 CONTINUE
564:       WORK( 1 ) = MAXWRK
565:       IWORK( 1 ) = LIWORK
566:       RWORK( 1 ) = LRWORK
567:       RETURN
568: *
569: *     End of ZGELSD
570: *
571:       END
572: