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] 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: *> \date December 2016
198: *
199: *> \ingroup doubleGEsolve
200: *
201: *> \par Contributors:
202: * ==================
203: *>
204: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
205: *> California at Berkeley, USA \n
206: *> Osni Marques, LBNL/NERSC, USA \n
207: *
208: * =====================================================================
209: SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
210: $ WORK, LWORK, IWORK, INFO )
211: *
212: * -- LAPACK driver routine (version 3.7.0) --
213: * -- LAPACK is a software package provided by Univ. of Tennessee, --
214: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215: * December 2016
216: *
217: * .. Scalar Arguments ..
218: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
219: DOUBLE PRECISION RCOND
220: * ..
221: * .. Array Arguments ..
222: INTEGER IWORK( * )
223: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
224: * ..
225: *
226: * =====================================================================
227: *
228: * .. Parameters ..
229: DOUBLE PRECISION ZERO, ONE, TWO
230: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
231: * ..
232: * .. Local Scalars ..
233: LOGICAL LQUERY
234: INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
235: $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
236: $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
237: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
238: * ..
239: * .. External Subroutines ..
240: EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD,
241: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA
242: * ..
243: * .. External Functions ..
244: INTEGER ILAENV
245: DOUBLE PRECISION DLAMCH, DLANGE
246: EXTERNAL ILAENV, DLAMCH, DLANGE
247: * ..
248: * .. Intrinsic Functions ..
249: INTRINSIC DBLE, INT, LOG, MAX, MIN
250: * ..
251: * .. Executable Statements ..
252: *
253: * Test the input arguments.
254: *
255: INFO = 0
256: MINMN = MIN( M, N )
257: MAXMN = MAX( M, N )
258: MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 )
259: LQUERY = ( LWORK.EQ.-1 )
260: IF( M.LT.0 ) THEN
261: INFO = -1
262: ELSE IF( N.LT.0 ) THEN
263: INFO = -2
264: ELSE IF( NRHS.LT.0 ) THEN
265: INFO = -3
266: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
267: INFO = -5
268: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
269: INFO = -7
270: END IF
271: *
272: SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
273: *
274: * Compute workspace.
275: * (Note: Comments in the code beginning "Workspace:" describe the
276: * minimal amount of workspace needed at that point in the code,
277: * as well as the preferred amount for good performance.
278: * NB refers to the optimal block size for the immediately
279: * following subroutine, as returned by ILAENV.)
280: *
281: MINWRK = 1
282: LIWORK = 1
283: MINMN = MAX( 1, MINMN )
284: NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) /
285: $ LOG( TWO ) ) + 1, 0 )
286: *
287: IF( INFO.EQ.0 ) THEN
288: MAXWRK = 0
289: LIWORK = 3*MINMN*NLVL + 11*MINMN
290: MM = M
291: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
292: *
293: * Path 1a - overdetermined, with many more rows than columns.
294: *
295: MM = N
296: MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N,
297: $ -1, -1 ) )
298: MAXWRK = MAX( MAXWRK, N+NRHS*
299: $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) )
300: END IF
301: IF( M.GE.N ) THEN
302: *
303: * Path 1 - overdetermined or exactly determined.
304: *
305: MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
306: $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) )
307: MAXWRK = MAX( MAXWRK, 3*N+NRHS*
308: $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) )
309: MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
310: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) )
311: WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
312: MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
313: MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
314: END IF
315: IF( N.GT.M ) THEN
316: WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
317: IF( N.GE.MNTHR ) THEN
318: *
319: * Path 2a - underdetermined, with many more columns
320: * than rows.
321: *
322: MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
323: MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
324: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
325: MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
326: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
327: MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
328: $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) )
329: IF( NRHS.GT.1 ) THEN
330: MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
331: ELSE
332: MAXWRK = MAX( MAXWRK, M*M+2*M )
333: END IF
334: MAXWRK = MAX( MAXWRK, M+NRHS*
335: $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) )
336: MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
337: ! XXX: Ensure the Path 2a case below is triggered. The workspace
338: ! calculation should use queries for all routines eventually.
339: MAXWRK = MAX( MAXWRK,
340: $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
341: ELSE
342: *
343: * Path 2 - remaining underdetermined cases.
344: *
345: MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N,
346: $ -1, -1 )
347: MAXWRK = MAX( MAXWRK, 3*M+NRHS*
348: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) )
349: MAXWRK = MAX( MAXWRK, 3*M+M*
350: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) )
351: MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
352: END IF
353: MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
354: END IF
355: MINWRK = MIN( MINWRK, MAXWRK )
356: WORK( 1 ) = MAXWRK
357: IWORK( 1 ) = LIWORK
358:
359: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
360: INFO = -12
361: END IF
362: END IF
363: *
364: IF( INFO.NE.0 ) THEN
365: CALL XERBLA( 'DGELSD', -INFO )
366: RETURN
367: ELSE IF( LQUERY ) THEN
368: GO TO 10
369: END IF
370: *
371: * Quick return if possible.
372: *
373: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
374: RANK = 0
375: RETURN
376: END IF
377: *
378: * Get machine parameters.
379: *
380: EPS = DLAMCH( 'P' )
381: SFMIN = DLAMCH( 'S' )
382: SMLNUM = SFMIN / EPS
383: BIGNUM = ONE / SMLNUM
384: CALL DLABAD( SMLNUM, BIGNUM )
385: *
386: * Scale A if max entry outside range [SMLNUM,BIGNUM].
387: *
388: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
389: IASCL = 0
390: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
391: *
392: * Scale matrix norm up to SMLNUM.
393: *
394: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
395: IASCL = 1
396: ELSE IF( ANRM.GT.BIGNUM ) THEN
397: *
398: * Scale matrix norm down to BIGNUM.
399: *
400: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
401: IASCL = 2
402: ELSE IF( ANRM.EQ.ZERO ) THEN
403: *
404: * Matrix all zero. Return zero solution.
405: *
406: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
407: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
408: RANK = 0
409: GO TO 10
410: END IF
411: *
412: * Scale B if max entry outside range [SMLNUM,BIGNUM].
413: *
414: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
415: IBSCL = 0
416: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
417: *
418: * Scale matrix norm up to SMLNUM.
419: *
420: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
421: IBSCL = 1
422: ELSE IF( BNRM.GT.BIGNUM ) THEN
423: *
424: * Scale matrix norm down to BIGNUM.
425: *
426: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
427: IBSCL = 2
428: END IF
429: *
430: * If M < N make sure certain entries of B are zero.
431: *
432: IF( M.LT.N )
433: $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
434: *
435: * Overdetermined case.
436: *
437: IF( M.GE.N ) THEN
438: *
439: * Path 1 - overdetermined or exactly determined.
440: *
441: MM = M
442: IF( M.GE.MNTHR ) THEN
443: *
444: * Path 1a - overdetermined, with many more rows than columns.
445: *
446: MM = N
447: ITAU = 1
448: NWORK = ITAU + N
449: *
450: * Compute A=Q*R.
451: * (Workspace: need 2*N, prefer N+N*NB)
452: *
453: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
454: $ LWORK-NWORK+1, INFO )
455: *
456: * Multiply B by transpose(Q).
457: * (Workspace: need N+NRHS, prefer N+NRHS*NB)
458: *
459: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
460: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
461: *
462: * Zero out below R.
463: *
464: IF( N.GT.1 ) THEN
465: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
466: END IF
467: END IF
468: *
469: IE = 1
470: ITAUQ = IE + N
471: ITAUP = ITAUQ + N
472: NWORK = ITAUP + N
473: *
474: * Bidiagonalize R in A.
475: * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
476: *
477: CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
478: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
479: $ INFO )
480: *
481: * Multiply B by transpose of left bidiagonalizing vectors of R.
482: * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
483: *
484: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
485: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
486: *
487: * Solve the bidiagonal least squares problem.
488: *
489: CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
490: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
491: IF( INFO.NE.0 ) THEN
492: GO TO 10
493: END IF
494: *
495: * Multiply B by right bidiagonalizing vectors of R.
496: *
497: CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
498: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
499: *
500: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
501: $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN
502: *
503: * Path 2a - underdetermined, with many more columns than rows
504: * and sufficient workspace for an efficient algorithm.
505: *
506: LDWORK = M
507: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
508: $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA
509: ITAU = 1
510: NWORK = M + 1
511: *
512: * Compute A=L*Q.
513: * (Workspace: need 2*M, prefer M+M*NB)
514: *
515: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
516: $ LWORK-NWORK+1, INFO )
517: IL = NWORK
518: *
519: * Copy L to WORK(IL), zeroing out above its diagonal.
520: *
521: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
522: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
523: $ LDWORK )
524: IE = IL + LDWORK*M
525: ITAUQ = IE + M
526: ITAUP = ITAUQ + M
527: NWORK = ITAUP + M
528: *
529: * Bidiagonalize L in WORK(IL).
530: * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
531: *
532: CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
533: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
534: $ LWORK-NWORK+1, INFO )
535: *
536: * Multiply B by transpose of left bidiagonalizing vectors of L.
537: * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
538: *
539: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
540: $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
541: $ LWORK-NWORK+1, INFO )
542: *
543: * Solve the bidiagonal least squares problem.
544: *
545: CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
546: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
547: IF( INFO.NE.0 ) THEN
548: GO TO 10
549: END IF
550: *
551: * Multiply B by right bidiagonalizing vectors of L.
552: *
553: CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
554: $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
555: $ LWORK-NWORK+1, INFO )
556: *
557: * Zero out below first M rows of B.
558: *
559: CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
560: NWORK = ITAU + M
561: *
562: * Multiply transpose(Q) by B.
563: * (Workspace: need M+NRHS, prefer M+NRHS*NB)
564: *
565: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
566: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
567: *
568: ELSE
569: *
570: * Path 2 - remaining underdetermined cases.
571: *
572: IE = 1
573: ITAUQ = IE + M
574: ITAUP = ITAUQ + M
575: NWORK = ITAUP + M
576: *
577: * Bidiagonalize A.
578: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
579: *
580: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
581: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
582: $ INFO )
583: *
584: * Multiply B by transpose of left bidiagonalizing vectors.
585: * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
586: *
587: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
588: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
589: *
590: * Solve the bidiagonal least squares problem.
591: *
592: CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
593: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
594: IF( INFO.NE.0 ) THEN
595: GO TO 10
596: END IF
597: *
598: * Multiply B by right bidiagonalizing vectors of A.
599: *
600: CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
601: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
602: *
603: END IF
604: *
605: * Undo scaling.
606: *
607: IF( IASCL.EQ.1 ) THEN
608: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
609: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
610: $ INFO )
611: ELSE IF( IASCL.EQ.2 ) THEN
612: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
613: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
614: $ INFO )
615: END IF
616: IF( IBSCL.EQ.1 ) THEN
617: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
618: ELSE IF( IBSCL.EQ.2 ) THEN
619: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
620: END IF
621: *
622: 10 CONTINUE
623: WORK( 1 ) = MAXWRK
624: IWORK( 1 ) = LIWORK
625: RETURN
626: *
627: * End of DGELSD
628: *
629: 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