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