1: *> \brief <b> DGELSS solves overdetermined or underdetermined systems 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 DGELSS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelss.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelss.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelss.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
22: * WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26: * DOUBLE PRECISION RCOND
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DGELSS computes the minimum norm solution to a real linear least
39: *> squares problem:
40: *>
41: *> Minimize 2-norm(| b - A*x |).
42: *>
43: *> using the singular value decomposition (SVD) of A. A is an M-by-N
44: *> matrix which may be rank-deficient.
45: *>
46: *> Several right hand side vectors b and solution vectors x can be
47: *> handled in a single call; they are stored as the columns of the
48: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
49: *> X.
50: *>
51: *> The effective rank of A is determined by treating as zero those
52: *> singular values which are less than RCOND times the largest singular
53: *> value.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] M
60: *> \verbatim
61: *> M is INTEGER
62: *> The number of rows of the matrix A. M >= 0.
63: *> \endverbatim
64: *>
65: *> \param[in] N
66: *> \verbatim
67: *> N is INTEGER
68: *> The number of columns of the matrix A. N >= 0.
69: *> \endverbatim
70: *>
71: *> \param[in] NRHS
72: *> \verbatim
73: *> NRHS is INTEGER
74: *> The number of right hand sides, i.e., the number of columns
75: *> of the matrices B and X. NRHS >= 0.
76: *> \endverbatim
77: *>
78: *> \param[in,out] A
79: *> \verbatim
80: *> A is DOUBLE PRECISION array, dimension (LDA,N)
81: *> On entry, the M-by-N matrix A.
82: *> On exit, the first min(m,n) rows of A are overwritten with
83: *> its right singular vectors, stored rowwise.
84: *> \endverbatim
85: *>
86: *> \param[in] LDA
87: *> \verbatim
88: *> LDA is INTEGER
89: *> The leading dimension of the array A. LDA >= max(1,M).
90: *> \endverbatim
91: *>
92: *> \param[in,out] B
93: *> \verbatim
94: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
95: *> On entry, the M-by-NRHS right hand side matrix B.
96: *> On exit, B is overwritten by the N-by-NRHS solution
97: *> matrix X. If m >= n and RANK = n, the residual
98: *> sum-of-squares for the solution in the i-th column is given
99: *> by the sum of squares of elements n+1:m in that column.
100: *> \endverbatim
101: *>
102: *> \param[in] LDB
103: *> \verbatim
104: *> LDB is INTEGER
105: *> The leading dimension of the array B. LDB >= max(1,max(M,N)).
106: *> \endverbatim
107: *>
108: *> \param[out] S
109: *> \verbatim
110: *> S is DOUBLE PRECISION array, dimension (min(M,N))
111: *> The singular values of A in decreasing order.
112: *> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
113: *> \endverbatim
114: *>
115: *> \param[in] RCOND
116: *> \verbatim
117: *> RCOND is DOUBLE PRECISION
118: *> RCOND is used to determine the effective rank of A.
119: *> Singular values S(i) <= RCOND*S(1) are treated as zero.
120: *> If RCOND < 0, machine precision is used instead.
121: *> \endverbatim
122: *>
123: *> \param[out] RANK
124: *> \verbatim
125: *> RANK is INTEGER
126: *> The effective rank of A, i.e., the number of singular values
127: *> which are greater than RCOND*S(1).
128: *> \endverbatim
129: *>
130: *> \param[out] WORK
131: *> \verbatim
132: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
133: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
134: *> \endverbatim
135: *>
136: *> \param[in] LWORK
137: *> \verbatim
138: *> LWORK is INTEGER
139: *> The dimension of the array WORK. LWORK >= 1, and also:
140: *> LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
141: *> For good performance, LWORK should generally be larger.
142: *>
143: *> If LWORK = -1, then a workspace query is assumed; the routine
144: *> only calculates the optimal size of the WORK array, returns
145: *> this value as the first entry of the WORK array, and no error
146: *> message related to LWORK is issued by XERBLA.
147: *> \endverbatim
148: *>
149: *> \param[out] INFO
150: *> \verbatim
151: *> INFO is INTEGER
152: *> = 0: successful exit
153: *> < 0: if INFO = -i, the i-th argument had an illegal value.
154: *> > 0: the algorithm for computing the SVD failed to converge;
155: *> if INFO = i, i off-diagonal elements of an intermediate
156: *> bidiagonal form did not converge to zero.
157: *> \endverbatim
158: *
159: * Authors:
160: * ========
161: *
162: *> \author Univ. of Tennessee
163: *> \author Univ. of California Berkeley
164: *> \author Univ. of Colorado Denver
165: *> \author NAG Ltd.
166: *
167: *> \ingroup doubleGEsolve
168: *
169: * =====================================================================
170: SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
171: $ WORK, LWORK, INFO )
172: *
173: * -- LAPACK driver routine --
174: * -- LAPACK is a software package provided by Univ. of Tennessee, --
175: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176: *
177: * .. Scalar Arguments ..
178: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
179: DOUBLE PRECISION RCOND
180: * ..
181: * .. Array Arguments ..
182: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
183: * ..
184: *
185: * =====================================================================
186: *
187: * .. Parameters ..
188: DOUBLE PRECISION ZERO, ONE
189: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
190: * ..
191: * .. Local Scalars ..
192: LOGICAL LQUERY
193: INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
194: $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
195: $ MAXWRK, MINMN, MINWRK, MM, MNTHR
196: INTEGER LWORK_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD,
197: $ LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ,
198: $ LWORK_DGELQF
199: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
200: * ..
201: * .. Local Arrays ..
202: DOUBLE PRECISION DUM( 1 )
203: * ..
204: * .. External Subroutines ..
205: EXTERNAL DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV,
206: $ DGEQRF, DLABAD, DLACPY, DLASCL, DLASET, DORGBR,
207: $ DORMBR, DORMLQ, DORMQR, DRSCL, XERBLA
208: * ..
209: * .. External Functions ..
210: INTEGER ILAENV
211: DOUBLE PRECISION DLAMCH, DLANGE
212: EXTERNAL ILAENV, DLAMCH, DLANGE
213: * ..
214: * .. Intrinsic Functions ..
215: INTRINSIC MAX, MIN
216: * ..
217: * .. Executable Statements ..
218: *
219: * Test the input arguments
220: *
221: INFO = 0
222: MINMN = MIN( M, N )
223: MAXMN = MAX( M, N )
224: LQUERY = ( LWORK.EQ.-1 )
225: IF( M.LT.0 ) THEN
226: INFO = -1
227: ELSE IF( N.LT.0 ) THEN
228: INFO = -2
229: ELSE IF( NRHS.LT.0 ) THEN
230: INFO = -3
231: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
232: INFO = -5
233: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
234: INFO = -7
235: END IF
236: *
237: * Compute workspace
238: * (Note: Comments in the code beginning "Workspace:" describe the
239: * minimal amount of workspace needed at that point in the code,
240: * as well as the preferred amount for good performance.
241: * NB refers to the optimal block size for the immediately
242: * following subroutine, as returned by ILAENV.)
243: *
244: IF( INFO.EQ.0 ) THEN
245: MINWRK = 1
246: MAXWRK = 1
247: IF( MINMN.GT.0 ) THEN
248: MM = M
249: MNTHR = ILAENV( 6, 'DGELSS', ' ', M, N, NRHS, -1 )
250: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
251: *
252: * Path 1a - overdetermined, with many more rows than
253: * columns
254: *
255: * Compute space needed for DGEQRF
256: CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
257: LWORK_DGEQRF = INT( DUM(1) )
258: * Compute space needed for DORMQR
259: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
260: $ LDB, DUM(1), -1, INFO )
261: LWORK_DORMQR = INT( DUM(1) )
262: MM = N
263: MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF )
264: MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR )
265: END IF
266: IF( M.GE.N ) THEN
267: *
268: * Path 1 - overdetermined or exactly determined
269: *
270: * Compute workspace needed for DBDSQR
271: *
272: BDSPAC = MAX( 1, 5*N )
273: * Compute space needed for DGEBRD
274: CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
275: $ DUM(1), DUM(1), -1, INFO )
276: LWORK_DGEBRD = INT( DUM(1) )
277: * Compute space needed for DORMBR
278: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
279: $ B, LDB, DUM(1), -1, INFO )
280: LWORK_DORMBR = INT( DUM(1) )
281: * Compute space needed for DORGBR
282: CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
283: $ DUM(1), -1, INFO )
284: LWORK_DORGBR = INT( DUM(1) )
285: * Compute total workspace needed
286: MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
287: MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR )
288: MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR )
289: MAXWRK = MAX( MAXWRK, BDSPAC )
290: MAXWRK = MAX( MAXWRK, N*NRHS )
291: MINWRK = MAX( 3*N + MM, 3*N + NRHS, BDSPAC )
292: MAXWRK = MAX( MINWRK, MAXWRK )
293: END IF
294: IF( N.GT.M ) THEN
295: *
296: * Compute workspace needed for DBDSQR
297: *
298: BDSPAC = MAX( 1, 5*M )
299: MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
300: IF( N.GE.MNTHR ) THEN
301: *
302: * Path 2a - underdetermined, with many more columns
303: * than rows
304: *
305: * Compute space needed for DGELQF
306: CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1),
307: $ -1, INFO )
308: LWORK_DGELQF = INT( DUM(1) )
309: * Compute space needed for DGEBRD
310: CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
311: $ DUM(1), DUM(1), -1, INFO )
312: LWORK_DGEBRD = INT( DUM(1) )
313: * Compute space needed for DORMBR
314: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA,
315: $ DUM(1), B, LDB, DUM(1), -1, INFO )
316: LWORK_DORMBR = INT( DUM(1) )
317: * Compute space needed for DORGBR
318: CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1),
319: $ DUM(1), -1, INFO )
320: LWORK_DORGBR = INT( DUM(1) )
321: * Compute space needed for DORMLQ
322: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
323: $ B, LDB, DUM(1), -1, INFO )
324: LWORK_DORMLQ = INT( DUM(1) )
325: * Compute total workspace needed
326: MAXWRK = M + LWORK_DGELQF
327: MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD )
328: MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORMBR )
329: MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DORGBR )
330: MAXWRK = MAX( MAXWRK, M*M + M + BDSPAC )
331: IF( NRHS.GT.1 ) THEN
332: MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
333: ELSE
334: MAXWRK = MAX( MAXWRK, M*M + 2*M )
335: END IF
336: MAXWRK = MAX( MAXWRK, M + LWORK_DORMLQ )
337: ELSE
338: *
339: * Path 2 - underdetermined
340: *
341: * Compute space needed for DGEBRD
342: CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
343: $ DUM(1), DUM(1), -1, INFO )
344: LWORK_DGEBRD = INT( DUM(1) )
345: * Compute space needed for DORMBR
346: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA,
347: $ DUM(1), B, LDB, DUM(1), -1, INFO )
348: LWORK_DORMBR = INT( DUM(1) )
349: * Compute space needed for DORGBR
350: CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1),
351: $ DUM(1), -1, INFO )
352: LWORK_DORGBR = INT( DUM(1) )
353: MAXWRK = 3*M + LWORK_DGEBRD
354: MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR )
355: MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR )
356: MAXWRK = MAX( MAXWRK, BDSPAC )
357: MAXWRK = MAX( MAXWRK, N*NRHS )
358: END IF
359: END IF
360: MAXWRK = MAX( MINWRK, MAXWRK )
361: END IF
362: WORK( 1 ) = MAXWRK
363: *
364: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
365: $ INFO = -12
366: END IF
367: *
368: IF( INFO.NE.0 ) THEN
369: CALL XERBLA( 'DGELSS', -INFO )
370: RETURN
371: ELSE IF( LQUERY ) THEN
372: RETURN
373: END IF
374: *
375: * Quick return if possible
376: *
377: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
378: RANK = 0
379: RETURN
380: END IF
381: *
382: * Get machine parameters
383: *
384: EPS = DLAMCH( 'P' )
385: SFMIN = DLAMCH( 'S' )
386: SMLNUM = SFMIN / EPS
387: BIGNUM = ONE / SMLNUM
388: CALL DLABAD( SMLNUM, BIGNUM )
389: *
390: * Scale A if max element outside range [SMLNUM,BIGNUM]
391: *
392: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
393: IASCL = 0
394: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
395: *
396: * Scale matrix norm up to SMLNUM
397: *
398: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
399: IASCL = 1
400: ELSE IF( ANRM.GT.BIGNUM ) THEN
401: *
402: * Scale matrix norm down to BIGNUM
403: *
404: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
405: IASCL = 2
406: ELSE IF( ANRM.EQ.ZERO ) THEN
407: *
408: * Matrix all zero. Return zero solution.
409: *
410: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
411: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
412: RANK = 0
413: GO TO 70
414: END IF
415: *
416: * Scale B if max element outside range [SMLNUM,BIGNUM]
417: *
418: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
419: IBSCL = 0
420: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
421: *
422: * Scale matrix norm up to SMLNUM
423: *
424: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
425: IBSCL = 1
426: ELSE IF( BNRM.GT.BIGNUM ) THEN
427: *
428: * Scale matrix norm down to BIGNUM
429: *
430: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
431: IBSCL = 2
432: END IF
433: *
434: * Overdetermined case
435: *
436: IF( M.GE.N ) THEN
437: *
438: * Path 1 - overdetermined or exactly determined
439: *
440: MM = M
441: IF( M.GE.MNTHR ) THEN
442: *
443: * Path 1a - overdetermined, with many more rows than columns
444: *
445: MM = N
446: ITAU = 1
447: IWORK = ITAU + N
448: *
449: * Compute A=Q*R
450: * (Workspace: need 2*N, prefer N+N*NB)
451: *
452: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
453: $ LWORK-IWORK+1, INFO )
454: *
455: * Multiply B by transpose(Q)
456: * (Workspace: need N+NRHS, prefer N+NRHS*NB)
457: *
458: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
459: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
460: *
461: * Zero out below R
462: *
463: IF( N.GT.1 )
464: $ CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
465: END IF
466: *
467: IE = 1
468: ITAUQ = IE + N
469: ITAUP = ITAUQ + N
470: IWORK = ITAUP + N
471: *
472: * Bidiagonalize R in A
473: * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
474: *
475: CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
476: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
477: $ INFO )
478: *
479: * Multiply B by transpose of left bidiagonalizing vectors of R
480: * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
481: *
482: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
483: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
484: *
485: * Generate right bidiagonalizing vectors of R in A
486: * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
487: *
488: CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
489: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
490: IWORK = IE + N
491: *
492: * Perform bidiagonal QR iteration
493: * multiply B by transpose of left singular vectors
494: * compute right singular vectors in A
495: * (Workspace: need BDSPAC)
496: *
497: CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
498: $ 1, B, LDB, WORK( IWORK ), INFO )
499: IF( INFO.NE.0 )
500: $ GO TO 70
501: *
502: * Multiply B by reciprocals of singular values
503: *
504: THR = MAX( RCOND*S( 1 ), SFMIN )
505: IF( RCOND.LT.ZERO )
506: $ THR = MAX( EPS*S( 1 ), SFMIN )
507: RANK = 0
508: DO 10 I = 1, N
509: IF( S( I ).GT.THR ) THEN
510: CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
511: RANK = RANK + 1
512: ELSE
513: CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
514: END IF
515: 10 CONTINUE
516: *
517: * Multiply B by right singular vectors
518: * (Workspace: need N, prefer N*NRHS)
519: *
520: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
521: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
522: $ WORK, LDB )
523: CALL DLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
524: ELSE IF( NRHS.GT.1 ) THEN
525: CHUNK = LWORK / N
526: DO 20 I = 1, NRHS, CHUNK
527: BL = MIN( NRHS-I+1, CHUNK )
528: CALL DGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B( 1, I ),
529: $ LDB, ZERO, WORK, N )
530: CALL DLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
531: 20 CONTINUE
532: ELSE
533: CALL DGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
534: CALL DCOPY( N, WORK, 1, B, 1 )
535: END IF
536: *
537: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
538: $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
539: *
540: * Path 2a - underdetermined, with many more columns than rows
541: * and sufficient workspace for an efficient algorithm
542: *
543: LDWORK = M
544: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
545: $ M*LDA+M+M*NRHS ) )LDWORK = LDA
546: ITAU = 1
547: IWORK = M + 1
548: *
549: * Compute A=L*Q
550: * (Workspace: need 2*M, prefer M+M*NB)
551: *
552: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
553: $ LWORK-IWORK+1, INFO )
554: IL = IWORK
555: *
556: * Copy L to WORK(IL), zeroing out above it
557: *
558: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
559: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
560: $ LDWORK )
561: IE = IL + LDWORK*M
562: ITAUQ = IE + M
563: ITAUP = ITAUQ + M
564: IWORK = ITAUP + M
565: *
566: * Bidiagonalize L in WORK(IL)
567: * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
568: *
569: CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
570: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
571: $ LWORK-IWORK+1, INFO )
572: *
573: * Multiply B by transpose of left bidiagonalizing vectors of L
574: * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
575: *
576: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
577: $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
578: $ LWORK-IWORK+1, INFO )
579: *
580: * Generate right bidiagonalizing vectors of R in WORK(IL)
581: * (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
582: *
583: CALL DORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
584: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
585: IWORK = IE + M
586: *
587: * Perform bidiagonal QR iteration,
588: * computing right singular vectors of L in WORK(IL) and
589: * multiplying B by transpose of left singular vectors
590: * (Workspace: need M*M+M+BDSPAC)
591: *
592: CALL DBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ),
593: $ LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO )
594: IF( INFO.NE.0 )
595: $ GO TO 70
596: *
597: * Multiply B by reciprocals of singular values
598: *
599: THR = MAX( RCOND*S( 1 ), SFMIN )
600: IF( RCOND.LT.ZERO )
601: $ THR = MAX( EPS*S( 1 ), SFMIN )
602: RANK = 0
603: DO 30 I = 1, M
604: IF( S( I ).GT.THR ) THEN
605: CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
606: RANK = RANK + 1
607: ELSE
608: CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
609: END IF
610: 30 CONTINUE
611: IWORK = IE
612: *
613: * Multiply B by right singular vectors of L in WORK(IL)
614: * (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
615: *
616: IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
617: CALL DGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK,
618: $ B, LDB, ZERO, WORK( IWORK ), LDB )
619: CALL DLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
620: ELSE IF( NRHS.GT.1 ) THEN
621: CHUNK = ( LWORK-IWORK+1 ) / M
622: DO 40 I = 1, NRHS, CHUNK
623: BL = MIN( NRHS-I+1, CHUNK )
624: CALL DGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK,
625: $ B( 1, I ), LDB, ZERO, WORK( IWORK ), M )
626: CALL DLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
627: $ LDB )
628: 40 CONTINUE
629: ELSE
630: CALL DGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ),
631: $ 1, ZERO, WORK( IWORK ), 1 )
632: CALL DCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
633: END IF
634: *
635: * Zero out below first M rows of B
636: *
637: CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
638: IWORK = ITAU + M
639: *
640: * Multiply transpose(Q) by B
641: * (Workspace: need M+NRHS, prefer M+NRHS*NB)
642: *
643: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
644: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
645: *
646: ELSE
647: *
648: * Path 2 - remaining underdetermined cases
649: *
650: IE = 1
651: ITAUQ = IE + M
652: ITAUP = ITAUQ + M
653: IWORK = ITAUP + M
654: *
655: * Bidiagonalize A
656: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
657: *
658: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
659: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
660: $ INFO )
661: *
662: * Multiply B by transpose of left bidiagonalizing vectors
663: * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
664: *
665: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
666: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
667: *
668: * Generate right bidiagonalizing vectors in A
669: * (Workspace: need 4*M, prefer 3*M+M*NB)
670: *
671: CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
672: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
673: IWORK = IE + M
674: *
675: * Perform bidiagonal QR iteration,
676: * computing right singular vectors of A in A and
677: * multiplying B by transpose of left singular vectors
678: * (Workspace: need BDSPAC)
679: *
680: CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, DUM,
681: $ 1, B, LDB, WORK( IWORK ), INFO )
682: IF( INFO.NE.0 )
683: $ GO TO 70
684: *
685: * Multiply B by reciprocals of singular values
686: *
687: THR = MAX( RCOND*S( 1 ), SFMIN )
688: IF( RCOND.LT.ZERO )
689: $ THR = MAX( EPS*S( 1 ), SFMIN )
690: RANK = 0
691: DO 50 I = 1, M
692: IF( S( I ).GT.THR ) THEN
693: CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
694: RANK = RANK + 1
695: ELSE
696: CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
697: END IF
698: 50 CONTINUE
699: *
700: * Multiply B by right singular vectors of A
701: * (Workspace: need N, prefer N*NRHS)
702: *
703: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
704: CALL DGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO,
705: $ WORK, LDB )
706: CALL DLACPY( 'F', N, NRHS, WORK, LDB, B, LDB )
707: ELSE IF( NRHS.GT.1 ) THEN
708: CHUNK = LWORK / N
709: DO 60 I = 1, NRHS, CHUNK
710: BL = MIN( NRHS-I+1, CHUNK )
711: CALL DGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ),
712: $ LDB, ZERO, WORK, N )
713: CALL DLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
714: 60 CONTINUE
715: ELSE
716: CALL DGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
717: CALL DCOPY( N, WORK, 1, B, 1 )
718: END IF
719: END IF
720: *
721: * Undo scaling
722: *
723: IF( IASCL.EQ.1 ) THEN
724: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
725: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
726: $ INFO )
727: ELSE IF( IASCL.EQ.2 ) THEN
728: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
729: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
730: $ INFO )
731: END IF
732: IF( IBSCL.EQ.1 ) THEN
733: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
734: ELSE IF( IBSCL.EQ.2 ) THEN
735: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
736: END IF
737: *
738: 70 CONTINUE
739: WORK( 1 ) = MAXWRK
740: RETURN
741: *
742: * End of DGELSS
743: *
744: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelss.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack