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