1 SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
2 $ WORK, LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
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 the modulus 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 (MAX(1,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 *
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 *
95 * RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
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
108 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
109 COMPLEX*16 CZERO, CONE
110 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
111 $ CONE = ( 1.0D+0, 0.0D+0 ) )
112 * ..
113 * .. Local Scalars ..
114 LOGICAL LQUERY
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 LQUERY = ( LWORK.EQ.-1 )
145 IF( M.LT.0 ) THEN
146 INFO = -1
147 ELSE IF( N.LT.0 ) THEN
148 INFO = -2
149 ELSE IF( NRHS.LT.0 ) THEN
150 INFO = -3
151 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
152 INFO = -5
153 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
154 INFO = -7
155 END IF
156 *
157 * Compute workspace
158 * (Note: Comments in the code beginning "Workspace:" describe the
159 * minimal amount of workspace needed at that point in the code,
160 * as well as the preferred amount for good performance.
161 * CWorkspace refers to complex workspace, and RWorkspace refers
162 * to real workspace. NB refers to the optimal block size for the
163 * immediately following subroutine, as returned by ILAENV.)
164 *
165 IF( INFO.EQ.0 ) THEN
166 MINWRK = 1
167 MAXWRK = 1
168 IF( MINMN.GT.0 ) THEN
169 MM = M
170 MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
171 IF( M.GE.N .AND. M.GE.MNTHR ) THEN
172 *
173 * Path 1a - overdetermined, with many more rows than
174 * columns
175 *
176 MM = N
177 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
178 $ N, -1, -1 ) )
179 MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
180 $ M, NRHS, N, -1 ) )
181 END IF
182 IF( M.GE.N ) THEN
183 *
184 * Path 1 - overdetermined or exactly determined
185 *
186 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
187 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
188 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
189 $ 'QLC', MM, NRHS, N, -1 ) )
190 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
191 $ 'ZUNGBR', 'P', N, N, N, -1 ) )
192 MAXWRK = MAX( MAXWRK, N*NRHS )
193 MINWRK = 2*N + MAX( NRHS, M )
194 END IF
195 IF( N.GT.M ) THEN
196 MINWRK = 2*M + MAX( NRHS, N )
197 IF( N.GE.MNTHR ) THEN
198 *
199 * Path 2a - underdetermined, with many more columns
200 * than rows
201 *
202 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
203 $ -1 )
204 MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
205 $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
206 MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
207 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
208 MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
209 $ 'ZUNGBR', 'P', M, M, M, -1 ) )
210 IF( NRHS.GT.1 ) THEN
211 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
212 ELSE
213 MAXWRK = MAX( MAXWRK, M*M + 2*M )
214 END IF
215 MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ',
216 $ 'LC', N, NRHS, M, -1 ) )
217 ELSE
218 *
219 * Path 2 - underdetermined
220 *
221 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
222 $ N, -1, -1 )
223 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
224 $ 'QLC', M, NRHS, M, -1 ) )
225 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR',
226 $ 'P', M, N, M, -1 ) )
227 MAXWRK = MAX( MAXWRK, N*NRHS )
228 END IF
229 END IF
230 MAXWRK = MAX( MINWRK, MAXWRK )
231 END IF
232 WORK( 1 ) = MAXWRK
233 *
234 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
235 $ INFO = -12
236 END IF
237 *
238 IF( INFO.NE.0 ) THEN
239 CALL XERBLA( 'ZGELSS', -INFO )
240 RETURN
241 ELSE IF( LQUERY ) THEN
242 RETURN
243 END IF
244 *
245 * Quick return if possible
246 *
247 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
248 RANK = 0
249 RETURN
250 END IF
251 *
252 * Get machine parameters
253 *
254 EPS = DLAMCH( 'P' )
255 SFMIN = DLAMCH( 'S' )
256 SMLNUM = SFMIN / EPS
257 BIGNUM = ONE / SMLNUM
258 CALL DLABAD( SMLNUM, BIGNUM )
259 *
260 * Scale A if max element outside range [SMLNUM,BIGNUM]
261 *
262 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
263 IASCL = 0
264 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
265 *
266 * Scale matrix norm up to SMLNUM
267 *
268 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
269 IASCL = 1
270 ELSE IF( ANRM.GT.BIGNUM ) THEN
271 *
272 * Scale matrix norm down to BIGNUM
273 *
274 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
275 IASCL = 2
276 ELSE IF( ANRM.EQ.ZERO ) THEN
277 *
278 * Matrix all zero. Return zero solution.
279 *
280 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
281 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
282 RANK = 0
283 GO TO 70
284 END IF
285 *
286 * Scale B if max element outside range [SMLNUM,BIGNUM]
287 *
288 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
289 IBSCL = 0
290 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
291 *
292 * Scale matrix norm up to SMLNUM
293 *
294 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
295 IBSCL = 1
296 ELSE IF( BNRM.GT.BIGNUM ) THEN
297 *
298 * Scale matrix norm down to BIGNUM
299 *
300 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
301 IBSCL = 2
302 END IF
303 *
304 * Overdetermined case
305 *
306 IF( M.GE.N ) THEN
307 *
308 * Path 1 - overdetermined or exactly determined
309 *
310 MM = M
311 IF( M.GE.MNTHR ) THEN
312 *
313 * Path 1a - overdetermined, with many more rows than columns
314 *
315 MM = N
316 ITAU = 1
317 IWORK = ITAU + N
318 *
319 * Compute A=Q*R
320 * (CWorkspace: need 2*N, prefer N+N*NB)
321 * (RWorkspace: none)
322 *
323 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
324 $ LWORK-IWORK+1, INFO )
325 *
326 * Multiply B by transpose(Q)
327 * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
328 * (RWorkspace: none)
329 *
330 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
331 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
332 *
333 * Zero out below R
334 *
335 IF( N.GT.1 )
336 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
337 $ LDA )
338 END IF
339 *
340 IE = 1
341 ITAUQ = 1
342 ITAUP = ITAUQ + N
343 IWORK = ITAUP + N
344 *
345 * Bidiagonalize R in A
346 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
347 * (RWorkspace: need N)
348 *
349 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
350 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
351 $ INFO )
352 *
353 * Multiply B by transpose of left bidiagonalizing vectors of R
354 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
355 * (RWorkspace: none)
356 *
357 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
358 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
359 *
360 * Generate right bidiagonalizing vectors of R in A
361 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
362 * (RWorkspace: none)
363 *
364 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
365 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
366 IRWORK = IE + N
367 *
368 * Perform bidiagonal QR iteration
369 * multiply B by transpose of left singular vectors
370 * compute right singular vectors in A
371 * (CWorkspace: none)
372 * (RWorkspace: need BDSPAC)
373 *
374 CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
375 $ 1, B, LDB, RWORK( IRWORK ), INFO )
376 IF( INFO.NE.0 )
377 $ GO TO 70
378 *
379 * Multiply B by reciprocals of singular values
380 *
381 THR = MAX( RCOND*S( 1 ), SFMIN )
382 IF( RCOND.LT.ZERO )
383 $ THR = MAX( EPS*S( 1 ), SFMIN )
384 RANK = 0
385 DO 10 I = 1, N
386 IF( S( I ).GT.THR ) THEN
387 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
388 RANK = RANK + 1
389 ELSE
390 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
391 END IF
392 10 CONTINUE
393 *
394 * Multiply B by right singular vectors
395 * (CWorkspace: need N, prefer N*NRHS)
396 * (RWorkspace: none)
397 *
398 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
399 CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
400 $ CZERO, WORK, LDB )
401 CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
402 ELSE IF( NRHS.GT.1 ) THEN
403 CHUNK = LWORK / N
404 DO 20 I = 1, NRHS, CHUNK
405 BL = MIN( NRHS-I+1, CHUNK )
406 CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
407 $ LDB, CZERO, WORK, N )
408 CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
409 20 CONTINUE
410 ELSE
411 CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
412 CALL ZCOPY( N, WORK, 1, B, 1 )
413 END IF
414 *
415 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
416 $ THEN
417 *
418 * Underdetermined case, M much less than N
419 *
420 * Path 2a - underdetermined, with many more columns than rows
421 * and sufficient workspace for an efficient algorithm
422 *
423 LDWORK = M
424 IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
425 $ LDWORK = LDA
426 ITAU = 1
427 IWORK = M + 1
428 *
429 * Compute A=L*Q
430 * (CWorkspace: need 2*M, prefer M+M*NB)
431 * (RWorkspace: none)
432 *
433 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
434 $ LWORK-IWORK+1, INFO )
435 IL = IWORK
436 *
437 * Copy L to WORK(IL), zeroing out above it
438 *
439 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
440 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
441 $ LDWORK )
442 IE = 1
443 ITAUQ = IL + LDWORK*M
444 ITAUP = ITAUQ + M
445 IWORK = ITAUP + M
446 *
447 * Bidiagonalize L in WORK(IL)
448 * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
449 * (RWorkspace: need M)
450 *
451 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
452 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
453 $ LWORK-IWORK+1, INFO )
454 *
455 * Multiply B by transpose of left bidiagonalizing vectors of L
456 * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
457 * (RWorkspace: none)
458 *
459 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
460 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
461 $ LWORK-IWORK+1, INFO )
462 *
463 * Generate right bidiagonalizing vectors of R in WORK(IL)
464 * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
465 * (RWorkspace: none)
466 *
467 CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
468 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
469 IRWORK = IE + M
470 *
471 * Perform bidiagonal QR iteration, computing right singular
472 * vectors of L in WORK(IL) and multiplying B by transpose of
473 * left singular vectors
474 * (CWorkspace: need M*M)
475 * (RWorkspace: need BDSPAC)
476 *
477 CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
478 $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
479 IF( INFO.NE.0 )
480 $ GO TO 70
481 *
482 * Multiply B by reciprocals of singular values
483 *
484 THR = MAX( RCOND*S( 1 ), SFMIN )
485 IF( RCOND.LT.ZERO )
486 $ THR = MAX( EPS*S( 1 ), SFMIN )
487 RANK = 0
488 DO 30 I = 1, M
489 IF( S( I ).GT.THR ) THEN
490 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
491 RANK = RANK + 1
492 ELSE
493 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
494 END IF
495 30 CONTINUE
496 IWORK = IL + M*LDWORK
497 *
498 * Multiply B by right singular vectors of L in WORK(IL)
499 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
500 * (RWorkspace: none)
501 *
502 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
503 CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
504 $ B, LDB, CZERO, WORK( IWORK ), LDB )
505 CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
506 ELSE IF( NRHS.GT.1 ) THEN
507 CHUNK = ( LWORK-IWORK+1 ) / M
508 DO 40 I = 1, NRHS, CHUNK
509 BL = MIN( NRHS-I+1, CHUNK )
510 CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
511 $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
512 CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
513 $ LDB )
514 40 CONTINUE
515 ELSE
516 CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
517 $ 1, CZERO, WORK( IWORK ), 1 )
518 CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
519 END IF
520 *
521 * Zero out below first M rows of B
522 *
523 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
524 IWORK = ITAU + M
525 *
526 * Multiply transpose(Q) by B
527 * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
528 * (RWorkspace: none)
529 *
530 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
531 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
532 *
533 ELSE
534 *
535 * Path 2 - remaining underdetermined cases
536 *
537 IE = 1
538 ITAUQ = 1
539 ITAUP = ITAUQ + M
540 IWORK = ITAUP + M
541 *
542 * Bidiagonalize A
543 * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
544 * (RWorkspace: need N)
545 *
546 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
547 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
548 $ INFO )
549 *
550 * Multiply B by transpose of left bidiagonalizing vectors
551 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
552 * (RWorkspace: none)
553 *
554 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
555 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
556 *
557 * Generate right bidiagonalizing vectors in A
558 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
559 * (RWorkspace: none)
560 *
561 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
562 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
563 IRWORK = IE + M
564 *
565 * Perform bidiagonal QR iteration,
566 * computing right singular vectors of A in A and
567 * multiplying B by transpose of left singular vectors
568 * (CWorkspace: none)
569 * (RWorkspace: need BDSPAC)
570 *
571 CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
572 $ 1, B, LDB, RWORK( IRWORK ), INFO )
573 IF( INFO.NE.0 )
574 $ GO TO 70
575 *
576 * Multiply B by reciprocals of singular values
577 *
578 THR = MAX( RCOND*S( 1 ), SFMIN )
579 IF( RCOND.LT.ZERO )
580 $ THR = MAX( EPS*S( 1 ), SFMIN )
581 RANK = 0
582 DO 50 I = 1, M
583 IF( S( I ).GT.THR ) THEN
584 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
585 RANK = RANK + 1
586 ELSE
587 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
588 END IF
589 50 CONTINUE
590 *
591 * Multiply B by right singular vectors of A
592 * (CWorkspace: need N, prefer N*NRHS)
593 * (RWorkspace: none)
594 *
595 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
596 CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
597 $ CZERO, WORK, LDB )
598 CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
599 ELSE IF( NRHS.GT.1 ) THEN
600 CHUNK = LWORK / N
601 DO 60 I = 1, NRHS, CHUNK
602 BL = MIN( NRHS-I+1, CHUNK )
603 CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
604 $ LDB, CZERO, WORK, N )
605 CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
606 60 CONTINUE
607 ELSE
608 CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
609 CALL ZCOPY( N, WORK, 1, B, 1 )
610 END IF
611 END IF
612 *
613 * Undo scaling
614 *
615 IF( IASCL.EQ.1 ) THEN
616 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
617 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
618 $ INFO )
619 ELSE IF( IASCL.EQ.2 ) THEN
620 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
621 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
622 $ INFO )
623 END IF
624 IF( IBSCL.EQ.1 ) THEN
625 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
626 ELSE IF( IBSCL.EQ.2 ) THEN
627 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
628 END IF
629 70 CONTINUE
630 WORK( 1 ) = MAXWRK
631 RETURN
632 *
633 * End of ZGELSS
634 *
635 END
2 $ WORK, LWORK, RWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
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 the modulus 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 (MAX(1,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 *
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 *
95 * RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
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
108 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
109 COMPLEX*16 CZERO, CONE
110 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
111 $ CONE = ( 1.0D+0, 0.0D+0 ) )
112 * ..
113 * .. Local Scalars ..
114 LOGICAL LQUERY
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 LQUERY = ( LWORK.EQ.-1 )
145 IF( M.LT.0 ) THEN
146 INFO = -1
147 ELSE IF( N.LT.0 ) THEN
148 INFO = -2
149 ELSE IF( NRHS.LT.0 ) THEN
150 INFO = -3
151 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
152 INFO = -5
153 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
154 INFO = -7
155 END IF
156 *
157 * Compute workspace
158 * (Note: Comments in the code beginning "Workspace:" describe the
159 * minimal amount of workspace needed at that point in the code,
160 * as well as the preferred amount for good performance.
161 * CWorkspace refers to complex workspace, and RWorkspace refers
162 * to real workspace. NB refers to the optimal block size for the
163 * immediately following subroutine, as returned by ILAENV.)
164 *
165 IF( INFO.EQ.0 ) THEN
166 MINWRK = 1
167 MAXWRK = 1
168 IF( MINMN.GT.0 ) THEN
169 MM = M
170 MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
171 IF( M.GE.N .AND. M.GE.MNTHR ) THEN
172 *
173 * Path 1a - overdetermined, with many more rows than
174 * columns
175 *
176 MM = N
177 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
178 $ N, -1, -1 ) )
179 MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
180 $ M, NRHS, N, -1 ) )
181 END IF
182 IF( M.GE.N ) THEN
183 *
184 * Path 1 - overdetermined or exactly determined
185 *
186 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
187 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
188 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
189 $ 'QLC', MM, NRHS, N, -1 ) )
190 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
191 $ 'ZUNGBR', 'P', N, N, N, -1 ) )
192 MAXWRK = MAX( MAXWRK, N*NRHS )
193 MINWRK = 2*N + MAX( NRHS, M )
194 END IF
195 IF( N.GT.M ) THEN
196 MINWRK = 2*M + MAX( NRHS, N )
197 IF( N.GE.MNTHR ) THEN
198 *
199 * Path 2a - underdetermined, with many more columns
200 * than rows
201 *
202 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
203 $ -1 )
204 MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
205 $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
206 MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
207 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
208 MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
209 $ 'ZUNGBR', 'P', M, M, M, -1 ) )
210 IF( NRHS.GT.1 ) THEN
211 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
212 ELSE
213 MAXWRK = MAX( MAXWRK, M*M + 2*M )
214 END IF
215 MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ',
216 $ 'LC', N, NRHS, M, -1 ) )
217 ELSE
218 *
219 * Path 2 - underdetermined
220 *
221 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
222 $ N, -1, -1 )
223 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
224 $ 'QLC', M, NRHS, M, -1 ) )
225 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR',
226 $ 'P', M, N, M, -1 ) )
227 MAXWRK = MAX( MAXWRK, N*NRHS )
228 END IF
229 END IF
230 MAXWRK = MAX( MINWRK, MAXWRK )
231 END IF
232 WORK( 1 ) = MAXWRK
233 *
234 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
235 $ INFO = -12
236 END IF
237 *
238 IF( INFO.NE.0 ) THEN
239 CALL XERBLA( 'ZGELSS', -INFO )
240 RETURN
241 ELSE IF( LQUERY ) THEN
242 RETURN
243 END IF
244 *
245 * Quick return if possible
246 *
247 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
248 RANK = 0
249 RETURN
250 END IF
251 *
252 * Get machine parameters
253 *
254 EPS = DLAMCH( 'P' )
255 SFMIN = DLAMCH( 'S' )
256 SMLNUM = SFMIN / EPS
257 BIGNUM = ONE / SMLNUM
258 CALL DLABAD( SMLNUM, BIGNUM )
259 *
260 * Scale A if max element outside range [SMLNUM,BIGNUM]
261 *
262 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
263 IASCL = 0
264 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
265 *
266 * Scale matrix norm up to SMLNUM
267 *
268 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
269 IASCL = 1
270 ELSE IF( ANRM.GT.BIGNUM ) THEN
271 *
272 * Scale matrix norm down to BIGNUM
273 *
274 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
275 IASCL = 2
276 ELSE IF( ANRM.EQ.ZERO ) THEN
277 *
278 * Matrix all zero. Return zero solution.
279 *
280 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
281 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
282 RANK = 0
283 GO TO 70
284 END IF
285 *
286 * Scale B if max element outside range [SMLNUM,BIGNUM]
287 *
288 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
289 IBSCL = 0
290 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
291 *
292 * Scale matrix norm up to SMLNUM
293 *
294 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
295 IBSCL = 1
296 ELSE IF( BNRM.GT.BIGNUM ) THEN
297 *
298 * Scale matrix norm down to BIGNUM
299 *
300 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
301 IBSCL = 2
302 END IF
303 *
304 * Overdetermined case
305 *
306 IF( M.GE.N ) THEN
307 *
308 * Path 1 - overdetermined or exactly determined
309 *
310 MM = M
311 IF( M.GE.MNTHR ) THEN
312 *
313 * Path 1a - overdetermined, with many more rows than columns
314 *
315 MM = N
316 ITAU = 1
317 IWORK = ITAU + N
318 *
319 * Compute A=Q*R
320 * (CWorkspace: need 2*N, prefer N+N*NB)
321 * (RWorkspace: none)
322 *
323 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
324 $ LWORK-IWORK+1, INFO )
325 *
326 * Multiply B by transpose(Q)
327 * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
328 * (RWorkspace: none)
329 *
330 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
331 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
332 *
333 * Zero out below R
334 *
335 IF( N.GT.1 )
336 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
337 $ LDA )
338 END IF
339 *
340 IE = 1
341 ITAUQ = 1
342 ITAUP = ITAUQ + N
343 IWORK = ITAUP + N
344 *
345 * Bidiagonalize R in A
346 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
347 * (RWorkspace: need N)
348 *
349 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
350 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
351 $ INFO )
352 *
353 * Multiply B by transpose of left bidiagonalizing vectors of R
354 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
355 * (RWorkspace: none)
356 *
357 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
358 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
359 *
360 * Generate right bidiagonalizing vectors of R in A
361 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
362 * (RWorkspace: none)
363 *
364 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
365 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
366 IRWORK = IE + N
367 *
368 * Perform bidiagonal QR iteration
369 * multiply B by transpose of left singular vectors
370 * compute right singular vectors in A
371 * (CWorkspace: none)
372 * (RWorkspace: need BDSPAC)
373 *
374 CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
375 $ 1, B, LDB, RWORK( IRWORK ), INFO )
376 IF( INFO.NE.0 )
377 $ GO TO 70
378 *
379 * Multiply B by reciprocals of singular values
380 *
381 THR = MAX( RCOND*S( 1 ), SFMIN )
382 IF( RCOND.LT.ZERO )
383 $ THR = MAX( EPS*S( 1 ), SFMIN )
384 RANK = 0
385 DO 10 I = 1, N
386 IF( S( I ).GT.THR ) THEN
387 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
388 RANK = RANK + 1
389 ELSE
390 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
391 END IF
392 10 CONTINUE
393 *
394 * Multiply B by right singular vectors
395 * (CWorkspace: need N, prefer N*NRHS)
396 * (RWorkspace: none)
397 *
398 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
399 CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
400 $ CZERO, WORK, LDB )
401 CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
402 ELSE IF( NRHS.GT.1 ) THEN
403 CHUNK = LWORK / N
404 DO 20 I = 1, NRHS, CHUNK
405 BL = MIN( NRHS-I+1, CHUNK )
406 CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
407 $ LDB, CZERO, WORK, N )
408 CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
409 20 CONTINUE
410 ELSE
411 CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
412 CALL ZCOPY( N, WORK, 1, B, 1 )
413 END IF
414 *
415 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
416 $ THEN
417 *
418 * Underdetermined case, M much less than N
419 *
420 * Path 2a - underdetermined, with many more columns than rows
421 * and sufficient workspace for an efficient algorithm
422 *
423 LDWORK = M
424 IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
425 $ LDWORK = LDA
426 ITAU = 1
427 IWORK = M + 1
428 *
429 * Compute A=L*Q
430 * (CWorkspace: need 2*M, prefer M+M*NB)
431 * (RWorkspace: none)
432 *
433 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
434 $ LWORK-IWORK+1, INFO )
435 IL = IWORK
436 *
437 * Copy L to WORK(IL), zeroing out above it
438 *
439 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
440 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
441 $ LDWORK )
442 IE = 1
443 ITAUQ = IL + LDWORK*M
444 ITAUP = ITAUQ + M
445 IWORK = ITAUP + M
446 *
447 * Bidiagonalize L in WORK(IL)
448 * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
449 * (RWorkspace: need M)
450 *
451 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
452 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
453 $ LWORK-IWORK+1, INFO )
454 *
455 * Multiply B by transpose of left bidiagonalizing vectors of L
456 * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
457 * (RWorkspace: none)
458 *
459 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
460 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
461 $ LWORK-IWORK+1, INFO )
462 *
463 * Generate right bidiagonalizing vectors of R in WORK(IL)
464 * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
465 * (RWorkspace: none)
466 *
467 CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
468 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
469 IRWORK = IE + M
470 *
471 * Perform bidiagonal QR iteration, computing right singular
472 * vectors of L in WORK(IL) and multiplying B by transpose of
473 * left singular vectors
474 * (CWorkspace: need M*M)
475 * (RWorkspace: need BDSPAC)
476 *
477 CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
478 $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
479 IF( INFO.NE.0 )
480 $ GO TO 70
481 *
482 * Multiply B by reciprocals of singular values
483 *
484 THR = MAX( RCOND*S( 1 ), SFMIN )
485 IF( RCOND.LT.ZERO )
486 $ THR = MAX( EPS*S( 1 ), SFMIN )
487 RANK = 0
488 DO 30 I = 1, M
489 IF( S( I ).GT.THR ) THEN
490 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
491 RANK = RANK + 1
492 ELSE
493 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
494 END IF
495 30 CONTINUE
496 IWORK = IL + M*LDWORK
497 *
498 * Multiply B by right singular vectors of L in WORK(IL)
499 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
500 * (RWorkspace: none)
501 *
502 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
503 CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
504 $ B, LDB, CZERO, WORK( IWORK ), LDB )
505 CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
506 ELSE IF( NRHS.GT.1 ) THEN
507 CHUNK = ( LWORK-IWORK+1 ) / M
508 DO 40 I = 1, NRHS, CHUNK
509 BL = MIN( NRHS-I+1, CHUNK )
510 CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
511 $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
512 CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
513 $ LDB )
514 40 CONTINUE
515 ELSE
516 CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
517 $ 1, CZERO, WORK( IWORK ), 1 )
518 CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
519 END IF
520 *
521 * Zero out below first M rows of B
522 *
523 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
524 IWORK = ITAU + M
525 *
526 * Multiply transpose(Q) by B
527 * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
528 * (RWorkspace: none)
529 *
530 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
531 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
532 *
533 ELSE
534 *
535 * Path 2 - remaining underdetermined cases
536 *
537 IE = 1
538 ITAUQ = 1
539 ITAUP = ITAUQ + M
540 IWORK = ITAUP + M
541 *
542 * Bidiagonalize A
543 * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
544 * (RWorkspace: need N)
545 *
546 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
547 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
548 $ INFO )
549 *
550 * Multiply B by transpose of left bidiagonalizing vectors
551 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
552 * (RWorkspace: none)
553 *
554 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
555 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
556 *
557 * Generate right bidiagonalizing vectors in A
558 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
559 * (RWorkspace: none)
560 *
561 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
562 $ WORK( IWORK ), LWORK-IWORK+1, INFO )
563 IRWORK = IE + M
564 *
565 * Perform bidiagonal QR iteration,
566 * computing right singular vectors of A in A and
567 * multiplying B by transpose of left singular vectors
568 * (CWorkspace: none)
569 * (RWorkspace: need BDSPAC)
570 *
571 CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
572 $ 1, B, LDB, RWORK( IRWORK ), INFO )
573 IF( INFO.NE.0 )
574 $ GO TO 70
575 *
576 * Multiply B by reciprocals of singular values
577 *
578 THR = MAX( RCOND*S( 1 ), SFMIN )
579 IF( RCOND.LT.ZERO )
580 $ THR = MAX( EPS*S( 1 ), SFMIN )
581 RANK = 0
582 DO 50 I = 1, M
583 IF( S( I ).GT.THR ) THEN
584 CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
585 RANK = RANK + 1
586 ELSE
587 CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
588 END IF
589 50 CONTINUE
590 *
591 * Multiply B by right singular vectors of A
592 * (CWorkspace: need N, prefer N*NRHS)
593 * (RWorkspace: none)
594 *
595 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
596 CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
597 $ CZERO, WORK, LDB )
598 CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
599 ELSE IF( NRHS.GT.1 ) THEN
600 CHUNK = LWORK / N
601 DO 60 I = 1, NRHS, CHUNK
602 BL = MIN( NRHS-I+1, CHUNK )
603 CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
604 $ LDB, CZERO, WORK, N )
605 CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
606 60 CONTINUE
607 ELSE
608 CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
609 CALL ZCOPY( N, WORK, 1, B, 1 )
610 END IF
611 END IF
612 *
613 * Undo scaling
614 *
615 IF( IASCL.EQ.1 ) THEN
616 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
617 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
618 $ INFO )
619 ELSE IF( IASCL.EQ.2 ) THEN
620 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
621 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
622 $ INFO )
623 END IF
624 IF( IBSCL.EQ.1 ) THEN
625 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
626 ELSE IF( IBSCL.EQ.2 ) THEN
627 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
628 END IF
629 70 CONTINUE
630 WORK( 1 ) = MAXWRK
631 RETURN
632 *
633 * End of ZGELSS
634 *
635 END