001:       SUBROUTINE CGELSD( 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:       REAL               RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               RWORK( * ), S( * )
016:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  CGELSD 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/output) COMPLEX 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 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) REAL 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) REAL
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 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) REAL 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:       REAL               ZERO, ONE, TWO
153:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
154:       COMPLEX            CZERO
155:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+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:       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
163: *     ..
164: *     .. External Subroutines ..
165:       EXTERNAL           CGEBRD, CGELQF, CGEQRF, CLACPY,
166:      $                   CLALSD, CLASCL, CLASET, CUNMBR,
167:      $                   CUNMLQ, CUNMQR, SLABAD, SLASCL,
168:      $                   SLASET, XERBLA
169: *     ..
170: *     .. External Functions ..
171:       INTEGER            ILAENV
172:       REAL               CLANGE, SLAMCH
173:       EXTERNAL           CLANGE, SLAMCH, ILAENV
174: *     ..
175: *     .. Intrinsic Functions ..
176:       INTRINSIC          INT, LOG, MAX, MIN, REAL
177: *     ..
178: *     .. Executable Statements ..
179: *
180: *     Test the input arguments.
181: *
182:       INFO = 0
183:       MINMN = MIN( M, N )
184:       MAXMN = MAX( M, N )
185:       LQUERY = ( LWORK.EQ.-1 )
186:       IF( M.LT.0 ) THEN
187:          INFO = -1
188:       ELSE IF( N.LT.0 ) THEN
189:          INFO = -2
190:       ELSE IF( NRHS.LT.0 ) THEN
191:          INFO = -3
192:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
193:          INFO = -5
194:       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
195:          INFO = -7
196:       END IF
197: *
198: *     Compute workspace.
199: *     (Note: Comments in the code beginning "Workspace:" describe the
200: *     minimal amount of workspace needed at that point in the code,
201: *     as well as the preferred amount for good performance.
202: *     NB refers to the optimal block size for the immediately
203: *     following subroutine, as returned by ILAENV.)
204: *
205:       IF( INFO.EQ.0 ) THEN
206:          MINWRK = 1
207:          MAXWRK = 1
208:          LIWORK = 1
209:          LRWORK = 1
210:          IF( MINMN.GT.0 ) THEN
211:             SMLSIZ = ILAENV( 9, 'CGELSD', ' ', 0, 0, 0, 0 )
212:             MNTHR = ILAENV( 6, 'CGELSD', ' ', M, N, NRHS, -1 )
213:             NLVL = MAX( INT( LOG( REAL( MINMN ) / REAL( SMLSIZ + 1 ) ) /
214:      $                  LOG( TWO ) ) + 1, 0 )
215:             LIWORK = 3*MINMN*NLVL + 11*MINMN
216:             MM = M
217:             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
218: *
219: *              Path 1a - overdetermined, with many more rows than
220: *                        columns.
221: *
222:                MM = N
223:                MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'CGEQRF', ' ', M, N,
224:      $                       -1, -1 ) )
225:                MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'CUNMQR', 'LC', M,
226:      $                       NRHS, N, -1 ) )
227:             END IF
228:             IF( M.GE.N ) THEN
229: *
230: *              Path 1 - overdetermined or exactly determined.
231: *
232:                LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
233:      $                  ( SMLSIZ + 1 )**2
234:                MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
235:      $                       'CGEBRD', ' ', MM, N, -1, -1 ) )
236:                MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR',
237:      $                       'QLC', MM, NRHS, N, -1 ) )
238:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
239:      $                       'CUNMBR', 'PLN', N, NRHS, N, -1 ) )
240:                MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
241:                MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
242:             END IF
243:             IF( N.GT.M ) THEN
244:                LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
245:      $                  ( SMLSIZ + 1 )**2
246:                IF( N.GE.MNTHR ) THEN
247: *
248: *                 Path 2a - underdetermined, with many more columns
249: *                           than rows.
250: *
251:                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
252:      $                     -1 )
253:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
254:      $                          'CGEBRD', ' ', M, M, -1, -1 ) )
255:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
256:      $                          'CUNMBR', 'QLC', M, NRHS, M, -1 ) )
257:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
258:      $                          'CUNMLQ', 'LC', N, NRHS, M, -1 ) )
259:                   IF( NRHS.GT.1 ) THEN
260:                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
261:                   ELSE
262:                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
263:                   END IF
264:                   MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
265: !     XXX: Ensure the Path 2a case below is triggered.  The workspace
266: !     calculation should use queries for all routines eventually.
267:                   MAXWRK = MAX( MAXWRK,
268:      $                 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
269:                ELSE
270: *
271: *                 Path 2 - underdetermined.
272: *
273:                   MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M,
274:      $                     N, -1, -1 )
275:                   MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR',
276:      $                          'QLC', M, NRHS, M, -1 ) )
277:                   MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNMBR',
278:      $                          'PLN', N, NRHS, M, -1 ) )
279:                   MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
280:                END IF
281:                MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
282:             END IF
283:          END IF
284:          MINWRK = MIN( MINWRK, MAXWRK )
285:          WORK( 1 ) = MAXWRK
286:          IWORK( 1 ) = LIWORK
287:          RWORK( 1 ) = LRWORK
288: *
289:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
290:             INFO = -12
291:          END IF
292:       END IF
293: *
294:       IF( INFO.NE.0 ) THEN
295:          CALL XERBLA( 'CGELSD', -INFO )
296:          RETURN
297:       ELSE IF( LQUERY ) THEN
298:          RETURN
299:       END IF
300: *
301: *     Quick return if possible.
302: *
303:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
304:          RANK = 0
305:          RETURN
306:       END IF
307: *
308: *     Get machine parameters.
309: *
310:       EPS = SLAMCH( 'P' )
311:       SFMIN = SLAMCH( 'S' )
312:       SMLNUM = SFMIN / EPS
313:       BIGNUM = ONE / SMLNUM
314:       CALL SLABAD( SMLNUM, BIGNUM )
315: *
316: *     Scale A if max entry outside range [SMLNUM,BIGNUM].
317: *
318:       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
319:       IASCL = 0
320:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
321: *
322: *        Scale matrix norm up to SMLNUM
323: *
324:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
325:          IASCL = 1
326:       ELSE IF( ANRM.GT.BIGNUM ) THEN
327: *
328: *        Scale matrix norm down to BIGNUM.
329: *
330:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
331:          IASCL = 2
332:       ELSE IF( ANRM.EQ.ZERO ) THEN
333: *
334: *        Matrix all zero. Return zero solution.
335: *
336:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
337:          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
338:          RANK = 0
339:          GO TO 10
340:       END IF
341: *
342: *     Scale B if max entry outside range [SMLNUM,BIGNUM].
343: *
344:       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
345:       IBSCL = 0
346:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
347: *
348: *        Scale matrix norm up to SMLNUM.
349: *
350:          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
351:          IBSCL = 1
352:       ELSE IF( BNRM.GT.BIGNUM ) THEN
353: *
354: *        Scale matrix norm down to BIGNUM.
355: *
356:          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
357:          IBSCL = 2
358:       END IF
359: *
360: *     If M < N make sure B(M+1:N,:) = 0
361: *
362:       IF( M.LT.N )
363:      $   CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
364: *
365: *     Overdetermined case.
366: *
367:       IF( M.GE.N ) THEN
368: *
369: *        Path 1 - overdetermined or exactly determined.
370: *
371:          MM = M
372:          IF( M.GE.MNTHR ) THEN
373: *
374: *           Path 1a - overdetermined, with many more rows than columns
375: *
376:             MM = N
377:             ITAU = 1
378:             NWORK = ITAU + N
379: *
380: *           Compute A=Q*R.
381: *           (RWorkspace: need N)
382: *           (CWorkspace: need N, prefer N*NB)
383: *
384:             CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
385:      $                   LWORK-NWORK+1, INFO )
386: *
387: *           Multiply B by transpose(Q).
388: *           (RWorkspace: need N)
389: *           (CWorkspace: need NRHS, prefer NRHS*NB)
390: *
391:             CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
392:      $                   LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
393: *
394: *           Zero out below R.
395: *
396:             IF( N.GT.1 ) THEN
397:                CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
398:      $                      LDA )
399:             END IF
400:          END IF
401: *
402:          ITAUQ = 1
403:          ITAUP = ITAUQ + N
404:          NWORK = ITAUP + N
405:          IE = 1
406:          NRWORK = IE + N
407: *
408: *        Bidiagonalize R in A.
409: *        (RWorkspace: need N)
410: *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
411: *
412:          CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
413:      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
414:      $                INFO )
415: *
416: *        Multiply B by transpose of left bidiagonalizing vectors of R.
417: *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
418: *
419:          CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
420:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
421: *
422: *        Solve the bidiagonal least squares problem.
423: *
424:          CALL CLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
425:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
426:      $                IWORK, INFO )
427:          IF( INFO.NE.0 ) THEN
428:             GO TO 10
429:          END IF
430: *
431: *        Multiply B by right bidiagonalizing vectors of R.
432: *
433:          CALL CUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
434:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
435: *
436:       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
437:      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
438: *
439: *        Path 2a - underdetermined, with many more columns than rows
440: *        and sufficient workspace for an efficient algorithm.
441: *
442:          LDWORK = M
443:          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
444:      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
445:          ITAU = 1
446:          NWORK = M + 1
447: *
448: *        Compute A=L*Q.
449: *        (CWorkspace: need 2*M, prefer M+M*NB)
450: *
451:          CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
452:      $                LWORK-NWORK+1, INFO )
453:          IL = NWORK
454: *
455: *        Copy L to WORK(IL), zeroing out above its diagonal.
456: *
457:          CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
458:          CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
459:      $                LDWORK )
460:          ITAUQ = IL + LDWORK*M
461:          ITAUP = ITAUQ + M
462:          NWORK = ITAUP + M
463:          IE = 1
464:          NRWORK = IE + M
465: *
466: *        Bidiagonalize L in WORK(IL).
467: *        (RWorkspace: need M)
468: *        (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
469: *
470:          CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
471:      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
472:      $                LWORK-NWORK+1, INFO )
473: *
474: *        Multiply B by transpose of left bidiagonalizing vectors of L.
475: *        (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
476: *
477:          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
478:      $                WORK( ITAUQ ), B, LDB, WORK( NWORK ),
479:      $                LWORK-NWORK+1, INFO )
480: *
481: *        Solve the bidiagonal least squares problem.
482: *
483:          CALL CLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
484:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
485:      $                IWORK, INFO )
486:          IF( INFO.NE.0 ) THEN
487:             GO TO 10
488:          END IF
489: *
490: *        Multiply B by right bidiagonalizing vectors of L.
491: *
492:          CALL CUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
493:      $                WORK( ITAUP ), B, LDB, WORK( NWORK ),
494:      $                LWORK-NWORK+1, INFO )
495: *
496: *        Zero out below first M rows of B.
497: *
498:          CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
499:          NWORK = ITAU + M
500: *
501: *        Multiply transpose(Q) by B.
502: *        (CWorkspace: need NRHS, prefer NRHS*NB)
503: *
504:          CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
505:      $                LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
506: *
507:       ELSE
508: *
509: *        Path 2 - remaining underdetermined cases.
510: *
511:          ITAUQ = 1
512:          ITAUP = ITAUQ + M
513:          NWORK = ITAUP + M
514:          IE = 1
515:          NRWORK = IE + M
516: *
517: *        Bidiagonalize A.
518: *        (RWorkspace: need M)
519: *        (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
520: *
521:          CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
522:      $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
523:      $                INFO )
524: *
525: *        Multiply B by transpose of left bidiagonalizing vectors.
526: *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
527: *
528:          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
529:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
530: *
531: *        Solve the bidiagonal least squares problem.
532: *
533:          CALL CLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
534:      $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
535:      $                IWORK, INFO )
536:          IF( INFO.NE.0 ) THEN
537:             GO TO 10
538:          END IF
539: *
540: *        Multiply B by right bidiagonalizing vectors of A.
541: *
542:          CALL CUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
543:      $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
544: *
545:       END IF
546: *
547: *     Undo scaling.
548: *
549:       IF( IASCL.EQ.1 ) THEN
550:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
551:          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
552:      $                INFO )
553:       ELSE IF( IASCL.EQ.2 ) THEN
554:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
555:          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
556:      $                INFO )
557:       END IF
558:       IF( IBSCL.EQ.1 ) THEN
559:          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
560:       ELSE IF( IBSCL.EQ.2 ) THEN
561:          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
562:       END IF
563: *
564:    10 CONTINUE
565:       WORK( 1 ) = MAXWRK
566:       IWORK( 1 ) = LIWORK
567:       RWORK( 1 ) = LRWORK
568:       RETURN
569: *
570: *     End of CGELSD
571: *
572:       END
573: