1: *> \brief <b> DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGELSD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
22: * WORK, LWORK, IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26: * DOUBLE PRECISION RCOND
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IWORK( * )
30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DGELSD computes the minimum-norm solution to a real linear least
40: *> squares problem:
41: *> minimize 2-norm(| b - A*x |)
42: *> using the singular value decomposition (SVD) of A. A is an M-by-N
43: *> matrix which may be rank-deficient.
44: *>
45: *> Several right hand side vectors b and solution vectors x can be
46: *> handled in a single call; they are stored as the columns of the
47: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48: *> matrix X.
49: *>
50: *> The problem is solved in three steps:
51: *> (1) Reduce the coefficient matrix A to bidiagonal form with
52: *> Householder transformations, reducing the original problem
53: *> into a "bidiagonal least squares problem" (BLS)
54: *> (2) Solve the BLS using a divide and conquer approach.
55: *> (3) Apply back all the Householder transformations to solve
56: *> the original least squares problem.
57: *>
58: *> The effective rank of A is determined by treating as zero those
59: *> singular values which are less than RCOND times the largest singular
60: *> value.
61: *>
62: *> The divide and conquer algorithm makes very mild assumptions about
63: *> floating point arithmetic. It will work on machines with a guard
64: *> digit in add/subtract, or on those binary machines without guard
65: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
66: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
67: *> without guard digits, but we know of none.
68: *> \endverbatim
69: *
70: * Arguments:
71: * ==========
72: *
73: *> \param[in] M
74: *> \verbatim
75: *> M is INTEGER
76: *> The number of rows of A. M >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in] N
80: *> \verbatim
81: *> N is INTEGER
82: *> The number of columns of A. N >= 0.
83: *> \endverbatim
84: *>
85: *> \param[in] NRHS
86: *> \verbatim
87: *> NRHS is INTEGER
88: *> The number of right hand sides, i.e., the number of columns
89: *> of the matrices B and X. NRHS >= 0.
90: *> \endverbatim
91: *>
92: *> \param[in,out] A
93: *> \verbatim
94: *> A is DOUBLE PRECISION array, dimension (LDA,N)
95: *> On entry, the M-by-N matrix A.
96: *> On exit, A has been destroyed.
97: *> \endverbatim
98: *>
99: *> \param[in] LDA
100: *> \verbatim
101: *> LDA is INTEGER
102: *> The leading dimension of the array A. LDA >= max(1,M).
103: *> \endverbatim
104: *>
105: *> \param[in,out] B
106: *> \verbatim
107: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
108: *> On entry, the M-by-NRHS right hand side matrix B.
109: *> On exit, B is overwritten by the N-by-NRHS solution
110: *> matrix X. If m >= n and RANK = n, the residual
111: *> sum-of-squares for the solution in the i-th column is given
112: *> by the sum of squares of elements n+1:m in that column.
113: *> \endverbatim
114: *>
115: *> \param[in] LDB
116: *> \verbatim
117: *> LDB is INTEGER
118: *> The leading dimension of the array B. LDB >= max(1,max(M,N)).
119: *> \endverbatim
120: *>
121: *> \param[out] S
122: *> \verbatim
123: *> S is DOUBLE PRECISION array, dimension (min(M,N))
124: *> The singular values of A in decreasing order.
125: *> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
126: *> \endverbatim
127: *>
128: *> \param[in] RCOND
129: *> \verbatim
130: *> RCOND is DOUBLE PRECISION
131: *> RCOND is used to determine the effective rank of A.
132: *> Singular values S(i) <= RCOND*S(1) are treated as zero.
133: *> If RCOND < 0, machine precision is used instead.
134: *> \endverbatim
135: *>
136: *> \param[out] RANK
137: *> \verbatim
138: *> RANK is INTEGER
139: *> The effective rank of A, i.e., the number of singular values
140: *> which are greater than RCOND*S(1).
141: *> \endverbatim
142: *>
143: *> \param[out] WORK
144: *> \verbatim
145: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
146: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
147: *> \endverbatim
148: *>
149: *> \param[in] LWORK
150: *> \verbatim
151: *> LWORK is INTEGER
152: *> The dimension of the array WORK. LWORK must be at least 1.
153: *> The exact minimum amount of workspace needed depends on M,
154: *> N and NRHS. As long as LWORK is at least
155: *> 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
156: *> if M is greater than or equal to N or
157: *> 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
158: *> if M is less than N, the code will execute correctly.
159: *> SMLSIZ is returned by ILAENV and is equal to the maximum
160: *> size of the subproblems at the bottom of the computation
161: *> tree (usually about 25), and
162: *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
163: *> For good performance, LWORK should generally be larger.
164: *>
165: *> If LWORK = -1, then a workspace query is assumed; the routine
166: *> only calculates the optimal size of the WORK array, returns
167: *> this value as the first entry of the WORK array, and no error
168: *> message related to LWORK is issued by XERBLA.
169: *> \endverbatim
170: *>
171: *> \param[out] IWORK
172: *> \verbatim
173: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
174: *> LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN),
175: *> where MINMN = MIN( M,N ).
176: *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
177: *> \endverbatim
178: *>
179: *> \param[out] INFO
180: *> \verbatim
181: *> INFO is INTEGER
182: *> = 0: successful exit
183: *> < 0: if INFO = -i, the i-th argument had an illegal value.
184: *> > 0: the algorithm for computing the SVD failed to converge;
185: *> if INFO = i, i off-diagonal elements of an intermediate
186: *> bidiagonal form did not converge to zero.
187: *> \endverbatim
188: *
189: * Authors:
190: * ========
191: *
192: *> \author Univ. of Tennessee
193: *> \author Univ. of California Berkeley
194: *> \author Univ. of Colorado Denver
195: *> \author NAG Ltd.
196: *
197: *> \ingroup doubleGEsolve
198: *
199: *> \par Contributors:
200: * ==================
201: *>
202: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
203: *> California at Berkeley, USA \n
204: *> Osni Marques, LBNL/NERSC, USA \n
205: *
206: * =====================================================================
207: SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
208: $ WORK, LWORK, IWORK, INFO )
209: *
210: * -- LAPACK driver routine --
211: * -- LAPACK is a software package provided by Univ. of Tennessee, --
212: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213: *
214: * .. Scalar Arguments ..
215: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
216: DOUBLE PRECISION RCOND
217: * ..
218: * .. Array Arguments ..
219: INTEGER IWORK( * )
220: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
221: * ..
222: *
223: * =====================================================================
224: *
225: * .. Parameters ..
226: DOUBLE PRECISION ZERO, ONE, TWO
227: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
228: * ..
229: * .. Local Scalars ..
230: LOGICAL LQUERY
231: INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
232: $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
233: $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
234: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
235: * ..
236: * .. External Subroutines ..
237: EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD,
238: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA
239: * ..
240: * .. External Functions ..
241: INTEGER ILAENV
242: DOUBLE PRECISION DLAMCH, DLANGE
243: EXTERNAL ILAENV, DLAMCH, DLANGE
244: * ..
245: * .. Intrinsic Functions ..
246: INTRINSIC DBLE, INT, LOG, MAX, MIN
247: * ..
248: * .. Executable Statements ..
249: *
250: * Test the input arguments.
251: *
252: INFO = 0
253: MINMN = MIN( M, N )
254: MAXMN = MAX( M, N )
255: MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 )
256: LQUERY = ( LWORK.EQ.-1 )
257: IF( M.LT.0 ) THEN
258: INFO = -1
259: ELSE IF( N.LT.0 ) THEN
260: INFO = -2
261: ELSE IF( NRHS.LT.0 ) THEN
262: INFO = -3
263: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
264: INFO = -5
265: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
266: INFO = -7
267: END IF
268: *
269: SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
270: *
271: * Compute workspace.
272: * (Note: Comments in the code beginning "Workspace:" describe the
273: * minimal amount of workspace needed at that point in the code,
274: * as well as the preferred amount for good performance.
275: * NB refers to the optimal block size for the immediately
276: * following subroutine, as returned by ILAENV.)
277: *
278: MINWRK = 1
279: LIWORK = 1
280: MINMN = MAX( 1, MINMN )
281: NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) /
282: $ LOG( TWO ) ) + 1, 0 )
283: *
284: IF( INFO.EQ.0 ) THEN
285: MAXWRK = 0
286: LIWORK = 3*MINMN*NLVL + 11*MINMN
287: MM = M
288: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
289: *
290: * Path 1a - overdetermined, with many more rows than columns.
291: *
292: MM = N
293: MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N,
294: $ -1, -1 ) )
295: MAXWRK = MAX( MAXWRK, N+NRHS*
296: $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) )
297: END IF
298: IF( M.GE.N ) THEN
299: *
300: * Path 1 - overdetermined or exactly determined.
301: *
302: MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
303: $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) )
304: MAXWRK = MAX( MAXWRK, 3*N+NRHS*
305: $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) )
306: MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
307: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) )
308: WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
309: MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
310: MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
311: END IF
312: IF( N.GT.M ) THEN
313: WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
314: IF( N.GE.MNTHR ) THEN
315: *
316: * Path 2a - underdetermined, with many more columns
317: * than rows.
318: *
319: MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
320: MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
321: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
322: MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
323: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
324: MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
325: $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) )
326: IF( NRHS.GT.1 ) THEN
327: MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
328: ELSE
329: MAXWRK = MAX( MAXWRK, M*M+2*M )
330: END IF
331: MAXWRK = MAX( MAXWRK, M+NRHS*
332: $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) )
333: MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
334: ! XXX: Ensure the Path 2a case below is triggered. The workspace
335: ! calculation should use queries for all routines eventually.
336: MAXWRK = MAX( MAXWRK,
337: $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
338: ELSE
339: *
340: * Path 2 - remaining underdetermined cases.
341: *
342: MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N,
343: $ -1, -1 )
344: MAXWRK = MAX( MAXWRK, 3*M+NRHS*
345: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) )
346: MAXWRK = MAX( MAXWRK, 3*M+M*
347: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) )
348: MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
349: END IF
350: MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
351: END IF
352: MINWRK = MIN( MINWRK, MAXWRK )
353: WORK( 1 ) = MAXWRK
354: IWORK( 1 ) = LIWORK
355:
356: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
357: INFO = -12
358: END IF
359: END IF
360: *
361: IF( INFO.NE.0 ) THEN
362: CALL XERBLA( 'DGELSD', -INFO )
363: RETURN
364: ELSE IF( LQUERY ) THEN
365: GO TO 10
366: END IF
367: *
368: * Quick return if possible.
369: *
370: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
371: RANK = 0
372: RETURN
373: END IF
374: *
375: * Get machine parameters.
376: *
377: EPS = DLAMCH( 'P' )
378: SFMIN = DLAMCH( 'S' )
379: SMLNUM = SFMIN / EPS
380: BIGNUM = ONE / SMLNUM
381: CALL DLABAD( SMLNUM, BIGNUM )
382: *
383: * Scale A if max entry outside range [SMLNUM,BIGNUM].
384: *
385: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
386: IASCL = 0
387: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
388: *
389: * Scale matrix norm up to SMLNUM.
390: *
391: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
392: IASCL = 1
393: ELSE IF( ANRM.GT.BIGNUM ) THEN
394: *
395: * Scale matrix norm down to BIGNUM.
396: *
397: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
398: IASCL = 2
399: ELSE IF( ANRM.EQ.ZERO ) THEN
400: *
401: * Matrix all zero. Return zero solution.
402: *
403: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
404: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
405: RANK = 0
406: GO TO 10
407: END IF
408: *
409: * Scale B if max entry outside range [SMLNUM,BIGNUM].
410: *
411: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
412: IBSCL = 0
413: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
414: *
415: * Scale matrix norm up to SMLNUM.
416: *
417: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
418: IBSCL = 1
419: ELSE IF( BNRM.GT.BIGNUM ) THEN
420: *
421: * Scale matrix norm down to BIGNUM.
422: *
423: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
424: IBSCL = 2
425: END IF
426: *
427: * If M < N make sure certain entries of B are zero.
428: *
429: IF( M.LT.N )
430: $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
431: *
432: * Overdetermined case.
433: *
434: IF( M.GE.N ) THEN
435: *
436: * Path 1 - overdetermined or exactly determined.
437: *
438: MM = M
439: IF( M.GE.MNTHR ) THEN
440: *
441: * Path 1a - overdetermined, with many more rows than columns.
442: *
443: MM = N
444: ITAU = 1
445: NWORK = ITAU + N
446: *
447: * Compute A=Q*R.
448: * (Workspace: need 2*N, prefer N+N*NB)
449: *
450: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
451: $ LWORK-NWORK+1, INFO )
452: *
453: * Multiply B by transpose(Q).
454: * (Workspace: need N+NRHS, prefer N+NRHS*NB)
455: *
456: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
457: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
458: *
459: * Zero out below R.
460: *
461: IF( N.GT.1 ) THEN
462: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
463: END IF
464: END IF
465: *
466: IE = 1
467: ITAUQ = IE + N
468: ITAUP = ITAUQ + N
469: NWORK = ITAUP + N
470: *
471: * Bidiagonalize R in A.
472: * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
473: *
474: CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
475: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
476: $ INFO )
477: *
478: * Multiply B by transpose of left bidiagonalizing vectors of R.
479: * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
480: *
481: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
482: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
483: *
484: * Solve the bidiagonal least squares problem.
485: *
486: CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
487: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
488: IF( INFO.NE.0 ) THEN
489: GO TO 10
490: END IF
491: *
492: * Multiply B by right bidiagonalizing vectors of R.
493: *
494: CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
495: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
496: *
497: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
498: $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN
499: *
500: * Path 2a - underdetermined, with many more columns than rows
501: * and sufficient workspace for an efficient algorithm.
502: *
503: LDWORK = M
504: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
505: $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA
506: ITAU = 1
507: NWORK = M + 1
508: *
509: * Compute A=L*Q.
510: * (Workspace: need 2*M, prefer M+M*NB)
511: *
512: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
513: $ LWORK-NWORK+1, INFO )
514: IL = NWORK
515: *
516: * Copy L to WORK(IL), zeroing out above its diagonal.
517: *
518: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
519: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
520: $ LDWORK )
521: IE = IL + LDWORK*M
522: ITAUQ = IE + M
523: ITAUP = ITAUQ + M
524: NWORK = ITAUP + M
525: *
526: * Bidiagonalize L in WORK(IL).
527: * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
528: *
529: CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
530: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
531: $ LWORK-NWORK+1, INFO )
532: *
533: * Multiply B by transpose of left bidiagonalizing vectors of L.
534: * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
535: *
536: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
537: $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
538: $ LWORK-NWORK+1, INFO )
539: *
540: * Solve the bidiagonal least squares problem.
541: *
542: CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
543: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
544: IF( INFO.NE.0 ) THEN
545: GO TO 10
546: END IF
547: *
548: * Multiply B by right bidiagonalizing vectors of L.
549: *
550: CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
551: $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
552: $ LWORK-NWORK+1, INFO )
553: *
554: * Zero out below first M rows of B.
555: *
556: CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
557: NWORK = ITAU + M
558: *
559: * Multiply transpose(Q) by B.
560: * (Workspace: need M+NRHS, prefer M+NRHS*NB)
561: *
562: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
563: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
564: *
565: ELSE
566: *
567: * Path 2 - remaining underdetermined cases.
568: *
569: IE = 1
570: ITAUQ = IE + M
571: ITAUP = ITAUQ + M
572: NWORK = ITAUP + M
573: *
574: * Bidiagonalize A.
575: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
576: *
577: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
578: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
579: $ INFO )
580: *
581: * Multiply B by transpose of left bidiagonalizing vectors.
582: * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
583: *
584: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
585: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
586: *
587: * Solve the bidiagonal least squares problem.
588: *
589: CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
590: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
591: IF( INFO.NE.0 ) THEN
592: GO TO 10
593: END IF
594: *
595: * Multiply B by right bidiagonalizing vectors of A.
596: *
597: CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
598: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
599: *
600: END IF
601: *
602: * Undo scaling.
603: *
604: IF( IASCL.EQ.1 ) THEN
605: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
606: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
607: $ INFO )
608: ELSE IF( IASCL.EQ.2 ) THEN
609: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
610: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
611: $ INFO )
612: END IF
613: IF( IBSCL.EQ.1 ) THEN
614: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
615: ELSE IF( IBSCL.EQ.2 ) THEN
616: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
617: END IF
618: *
619: 10 CONTINUE
620: WORK( 1 ) = MAXWRK
621: IWORK( 1 ) = LIWORK
622: RETURN
623: *
624: * End of DGELSD
625: *
626: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack