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