1 SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
2 $ WORK, LWORK, RWORK, IWORK, 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 INTEGER IWORK( * )
15 DOUBLE PRECISION RWORK( * ), S( * )
16 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZGELSD computes the minimum-norm solution to a real linear least
23 * squares problem:
24 * minimize 2-norm(| b - A*x |)
25 * using the singular value decomposition (SVD) of A. A is an M-by-N
26 * matrix which may be rank-deficient.
27 *
28 * Several right hand side vectors b and solution vectors x can be
29 * handled in a single call; they are stored as the columns of the
30 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
31 * matrix X.
32 *
33 * The problem is solved in three steps:
34 * (1) Reduce the coefficient matrix A to bidiagonal form with
35 * Householder tranformations, reducing the original problem
36 * into a "bidiagonal least squares problem" (BLS)
37 * (2) Solve the BLS using a divide and conquer approach.
38 * (3) Apply back all the Householder tranformations to solve
39 * the original least squares problem.
40 *
41 * The effective rank of A is determined by treating as zero those
42 * singular values which are less than RCOND times the largest singular
43 * value.
44 *
45 * The divide and conquer algorithm makes very mild assumptions about
46 * floating point arithmetic. It will work on machines with a guard
47 * digit in add/subtract, or on those binary machines without guard
48 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
49 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
50 * without guard digits, but we know of none.
51 *
52 * Arguments
53 * =========
54 *
55 * M (input) INTEGER
56 * The number of rows of the matrix A. M >= 0.
57 *
58 * N (input) INTEGER
59 * The number of columns of the matrix A. N >= 0.
60 *
61 * NRHS (input) INTEGER
62 * The number of right hand sides, i.e., the number of columns
63 * of the matrices B and X. NRHS >= 0.
64 *
65 * A (input) COMPLEX*16 array, dimension (LDA,N)
66 * On entry, the M-by-N matrix A.
67 * On exit, A has been destroyed.
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,M).
71 *
72 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
73 * On entry, the M-by-NRHS right hand side matrix B.
74 * On exit, B is overwritten by the N-by-NRHS solution matrix X.
75 * If m >= n and RANK = n, the residual sum-of-squares for
76 * the solution in the i-th column is given by the sum of
77 * squares of the modulus of elements n+1:m in that column.
78 *
79 * LDB (input) INTEGER
80 * The leading dimension of the array B. LDB >= max(1,M,N).
81 *
82 * S (output) DOUBLE PRECISION array, dimension (min(M,N))
83 * The singular values of A in decreasing order.
84 * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
85 *
86 * RCOND (input) DOUBLE PRECISION
87 * RCOND is used to determine the effective rank of A.
88 * Singular values S(i) <= RCOND*S(1) are treated as zero.
89 * If RCOND < 0, machine precision is used instead.
90 *
91 * RANK (output) INTEGER
92 * The effective rank of A, i.e., the number of singular values
93 * which are greater than RCOND*S(1).
94 *
95 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
96 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
97 *
98 * LWORK (input) INTEGER
99 * The dimension of the array WORK. LWORK must be at least 1.
100 * The exact minimum amount of workspace needed depends on M,
101 * N and NRHS. As long as LWORK is at least
102 * 2*N + N*NRHS
103 * if M is greater than or equal to N or
104 * 2*M + M*NRHS
105 * if M is less than N, the code will execute correctly.
106 * For good performance, LWORK should generally be larger.
107 *
108 * If LWORK = -1, then a workspace query is assumed; the routine
109 * only calculates the optimal size of the array WORK and the
110 * minimum sizes of the arrays RWORK and IWORK, and returns
111 * these values as the first entries of the WORK, RWORK and
112 * IWORK arrays, and no error message related to LWORK is issued
113 * by XERBLA.
114 *
115 * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
116 * LRWORK >=
117 * 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
118 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
119 * if M is greater than or equal to N or
120 * 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
121 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
122 * if M is less than N, the code will execute correctly.
123 * SMLSIZ is returned by ILAENV and is equal to the maximum
124 * size of the subproblems at the bottom of the computation
125 * tree (usually about 25), and
126 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
127 * On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
128 *
129 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
130 * LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
131 * where MINMN = MIN( M,N ).
132 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
133 *
134 * INFO (output) INTEGER
135 * = 0: successful exit
136 * < 0: if INFO = -i, the i-th argument had an illegal value.
137 * > 0: the algorithm for computing the SVD failed to converge;
138 * if INFO = i, i off-diagonal elements of an intermediate
139 * bidiagonal form did not converge to zero.
140 *
141 * Further Details
142 * ===============
143 *
144 * Based on contributions by
145 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
146 * California at Berkeley, USA
147 * Osni Marques, LBNL/NERSC, USA
148 *
149 * =====================================================================
150 *
151 * .. Parameters ..
152 DOUBLE PRECISION ZERO, ONE, TWO
153 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
154 COMPLEX*16 CZERO
155 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
156 * ..
157 * .. Local Scalars ..
158 LOGICAL LQUERY
159 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
160 $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
161 $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
162 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF,
166 $ ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR,
167 $ ZUNMLQ, ZUNMQR
168 * ..
169 * .. External Functions ..
170 INTEGER ILAENV
171 DOUBLE PRECISION DLAMCH, ZLANGE
172 EXTERNAL ILAENV, DLAMCH, ZLANGE
173 * ..
174 * .. Intrinsic Functions ..
175 INTRINSIC INT, LOG, MAX, MIN, DBLE
176 * ..
177 * .. Executable Statements ..
178 *
179 * Test the input arguments.
180 *
181 INFO = 0
182 MINMN = MIN( M, N )
183 MAXMN = MAX( M, N )
184 LQUERY = ( LWORK.EQ.-1 )
185 IF( M.LT.0 ) THEN
186 INFO = -1
187 ELSE IF( N.LT.0 ) THEN
188 INFO = -2
189 ELSE IF( NRHS.LT.0 ) THEN
190 INFO = -3
191 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
192 INFO = -5
193 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
194 INFO = -7
195 END IF
196 *
197 * Compute workspace.
198 * (Note: Comments in the code beginning "Workspace:" describe the
199 * minimal amount of workspace needed at that point in the code,
200 * as well as the preferred amount for good performance.
201 * NB refers to the optimal block size for the immediately
202 * following subroutine, as returned by ILAENV.)
203 *
204 IF( INFO.EQ.0 ) THEN
205 MINWRK = 1
206 MAXWRK = 1
207 LIWORK = 1
208 LRWORK = 1
209 IF( MINMN.GT.0 ) THEN
210 SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 )
211 MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 )
212 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) /
213 $ LOG( TWO ) ) + 1, 0 )
214 LIWORK = 3*MINMN*NLVL + 11*MINMN
215 MM = M
216 IF( M.GE.N .AND. M.GE.MNTHR ) THEN
217 *
218 * Path 1a - overdetermined, with many more rows than
219 * columns.
220 *
221 MM = N
222 MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N,
223 $ -1, -1 ) )
224 MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M,
225 $ NRHS, N, -1 ) )
226 END IF
227 IF( M.GE.N ) THEN
228 *
229 * Path 1 - overdetermined or exactly determined.
230 *
231 LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
232 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
233 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
234 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
235 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
236 $ 'QLC', MM, NRHS, N, -1 ) )
237 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
238 $ 'ZUNMBR', 'PLN', N, NRHS, N, -1 ) )
239 MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
240 MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
241 END IF
242 IF( N.GT.M ) THEN
243 LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
244 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
245 IF( N.GE.MNTHR ) THEN
246 *
247 * Path 2a - underdetermined, with many more columns
248 * than rows.
249 *
250 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
251 $ -1 )
252 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
253 $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
254 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
255 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
256 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
257 $ 'ZUNMLQ', 'LC', N, NRHS, M, -1 ) )
258 IF( NRHS.GT.1 ) THEN
259 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
260 ELSE
261 MAXWRK = MAX( MAXWRK, M*M + 2*M )
262 END IF
263 MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
264 ! XXX: Ensure the Path 2a case below is triggered. The workspace
265 ! calculation should use queries for all routines eventually.
266 MAXWRK = MAX( MAXWRK,
267 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
268 ELSE
269 *
270 * Path 2 - underdetermined.
271 *
272 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
273 $ N, -1, -1 )
274 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
275 $ 'QLC', M, NRHS, M, -1 ) )
276 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR',
277 $ 'PLN', N, NRHS, M, -1 ) )
278 MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
279 END IF
280 MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
281 END IF
282 END IF
283 MINWRK = MIN( MINWRK, MAXWRK )
284 WORK( 1 ) = MAXWRK
285 IWORK( 1 ) = LIWORK
286 RWORK( 1 ) = LRWORK
287 *
288 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
289 INFO = -12
290 END IF
291 END IF
292 *
293 IF( INFO.NE.0 ) THEN
294 CALL XERBLA( 'ZGELSD', -INFO )
295 RETURN
296 ELSE IF( LQUERY ) THEN
297 RETURN
298 END IF
299 *
300 * Quick return if possible.
301 *
302 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
303 RANK = 0
304 RETURN
305 END IF
306 *
307 * Get machine parameters.
308 *
309 EPS = DLAMCH( 'P' )
310 SFMIN = DLAMCH( 'S' )
311 SMLNUM = SFMIN / EPS
312 BIGNUM = ONE / SMLNUM
313 CALL DLABAD( SMLNUM, BIGNUM )
314 *
315 * Scale A if max entry outside range [SMLNUM,BIGNUM].
316 *
317 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
318 IASCL = 0
319 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
320 *
321 * Scale matrix norm up to SMLNUM
322 *
323 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
324 IASCL = 1
325 ELSE IF( ANRM.GT.BIGNUM ) THEN
326 *
327 * Scale matrix norm down to BIGNUM.
328 *
329 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
330 IASCL = 2
331 ELSE IF( ANRM.EQ.ZERO ) THEN
332 *
333 * Matrix all zero. Return zero solution.
334 *
335 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
336 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
337 RANK = 0
338 GO TO 10
339 END IF
340 *
341 * Scale B if max entry outside range [SMLNUM,BIGNUM].
342 *
343 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
344 IBSCL = 0
345 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
346 *
347 * Scale matrix norm up to SMLNUM.
348 *
349 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
350 IBSCL = 1
351 ELSE IF( BNRM.GT.BIGNUM ) THEN
352 *
353 * Scale matrix norm down to BIGNUM.
354 *
355 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
356 IBSCL = 2
357 END IF
358 *
359 * If M < N make sure B(M+1:N,:) = 0
360 *
361 IF( M.LT.N )
362 $ CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
363 *
364 * Overdetermined case.
365 *
366 IF( M.GE.N ) THEN
367 *
368 * Path 1 - overdetermined or exactly determined.
369 *
370 MM = M
371 IF( M.GE.MNTHR ) THEN
372 *
373 * Path 1a - overdetermined, with many more rows than columns
374 *
375 MM = N
376 ITAU = 1
377 NWORK = ITAU + N
378 *
379 * Compute A=Q*R.
380 * (RWorkspace: need N)
381 * (CWorkspace: need N, prefer N*NB)
382 *
383 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
384 $ LWORK-NWORK+1, INFO )
385 *
386 * Multiply B by transpose(Q).
387 * (RWorkspace: need N)
388 * (CWorkspace: need NRHS, prefer NRHS*NB)
389 *
390 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
391 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
392 *
393 * Zero out below R.
394 *
395 IF( N.GT.1 ) THEN
396 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
397 $ LDA )
398 END IF
399 END IF
400 *
401 ITAUQ = 1
402 ITAUP = ITAUQ + N
403 NWORK = ITAUP + N
404 IE = 1
405 NRWORK = IE + N
406 *
407 * Bidiagonalize R in A.
408 * (RWorkspace: need N)
409 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
410 *
411 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
412 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
413 $ INFO )
414 *
415 * Multiply B by transpose of left bidiagonalizing vectors of R.
416 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
417 *
418 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
419 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
420 *
421 * Solve the bidiagonal least squares problem.
422 *
423 CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
424 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
425 $ IWORK, INFO )
426 IF( INFO.NE.0 ) THEN
427 GO TO 10
428 END IF
429 *
430 * Multiply B by right bidiagonalizing vectors of R.
431 *
432 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
433 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
434 *
435 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
436 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
437 *
438 * Path 2a - underdetermined, with many more columns than rows
439 * and sufficient workspace for an efficient algorithm.
440 *
441 LDWORK = M
442 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
443 $ M*LDA+M+M*NRHS ) )LDWORK = LDA
444 ITAU = 1
445 NWORK = M + 1
446 *
447 * Compute A=L*Q.
448 * (CWorkspace: need 2*M, prefer M+M*NB)
449 *
450 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
451 $ LWORK-NWORK+1, INFO )
452 IL = NWORK
453 *
454 * Copy L to WORK(IL), zeroing out above its diagonal.
455 *
456 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
457 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
458 $ LDWORK )
459 ITAUQ = IL + LDWORK*M
460 ITAUP = ITAUQ + M
461 NWORK = ITAUP + M
462 IE = 1
463 NRWORK = IE + M
464 *
465 * Bidiagonalize L in WORK(IL).
466 * (RWorkspace: need M)
467 * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
468 *
469 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
470 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
471 $ LWORK-NWORK+1, INFO )
472 *
473 * Multiply B by transpose of left bidiagonalizing vectors of L.
474 * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
475 *
476 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
477 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
478 $ LWORK-NWORK+1, INFO )
479 *
480 * Solve the bidiagonal least squares problem.
481 *
482 CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
483 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
484 $ IWORK, INFO )
485 IF( INFO.NE.0 ) THEN
486 GO TO 10
487 END IF
488 *
489 * Multiply B by right bidiagonalizing vectors of L.
490 *
491 CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
492 $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
493 $ LWORK-NWORK+1, INFO )
494 *
495 * Zero out below first M rows of B.
496 *
497 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
498 NWORK = ITAU + M
499 *
500 * Multiply transpose(Q) by B.
501 * (CWorkspace: need NRHS, prefer NRHS*NB)
502 *
503 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
504 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
505 *
506 ELSE
507 *
508 * Path 2 - remaining underdetermined cases.
509 *
510 ITAUQ = 1
511 ITAUP = ITAUQ + M
512 NWORK = ITAUP + M
513 IE = 1
514 NRWORK = IE + M
515 *
516 * Bidiagonalize A.
517 * (RWorkspace: need M)
518 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
519 *
520 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
521 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
522 $ INFO )
523 *
524 * Multiply B by transpose of left bidiagonalizing vectors.
525 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
526 *
527 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
528 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
529 *
530 * Solve the bidiagonal least squares problem.
531 *
532 CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
533 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
534 $ IWORK, INFO )
535 IF( INFO.NE.0 ) THEN
536 GO TO 10
537 END IF
538 *
539 * Multiply B by right bidiagonalizing vectors of A.
540 *
541 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
542 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
543 *
544 END IF
545 *
546 * Undo scaling.
547 *
548 IF( IASCL.EQ.1 ) THEN
549 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
550 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
551 $ INFO )
552 ELSE IF( IASCL.EQ.2 ) THEN
553 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
554 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
555 $ INFO )
556 END IF
557 IF( IBSCL.EQ.1 ) THEN
558 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
559 ELSE IF( IBSCL.EQ.2 ) THEN
560 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
561 END IF
562 *
563 10 CONTINUE
564 WORK( 1 ) = MAXWRK
565 IWORK( 1 ) = LIWORK
566 RWORK( 1 ) = LRWORK
567 RETURN
568 *
569 * End of ZGELSD
570 *
571 END
2 $ WORK, LWORK, RWORK, IWORK, 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 INTEGER IWORK( * )
15 DOUBLE PRECISION RWORK( * ), S( * )
16 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZGELSD computes the minimum-norm solution to a real linear least
23 * squares problem:
24 * minimize 2-norm(| b - A*x |)
25 * using the singular value decomposition (SVD) of A. A is an M-by-N
26 * matrix which may be rank-deficient.
27 *
28 * Several right hand side vectors b and solution vectors x can be
29 * handled in a single call; they are stored as the columns of the
30 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
31 * matrix X.
32 *
33 * The problem is solved in three steps:
34 * (1) Reduce the coefficient matrix A to bidiagonal form with
35 * Householder tranformations, reducing the original problem
36 * into a "bidiagonal least squares problem" (BLS)
37 * (2) Solve the BLS using a divide and conquer approach.
38 * (3) Apply back all the Householder tranformations to solve
39 * the original least squares problem.
40 *
41 * The effective rank of A is determined by treating as zero those
42 * singular values which are less than RCOND times the largest singular
43 * value.
44 *
45 * The divide and conquer algorithm makes very mild assumptions about
46 * floating point arithmetic. It will work on machines with a guard
47 * digit in add/subtract, or on those binary machines without guard
48 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
49 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
50 * without guard digits, but we know of none.
51 *
52 * Arguments
53 * =========
54 *
55 * M (input) INTEGER
56 * The number of rows of the matrix A. M >= 0.
57 *
58 * N (input) INTEGER
59 * The number of columns of the matrix A. N >= 0.
60 *
61 * NRHS (input) INTEGER
62 * The number of right hand sides, i.e., the number of columns
63 * of the matrices B and X. NRHS >= 0.
64 *
65 * A (input) COMPLEX*16 array, dimension (LDA,N)
66 * On entry, the M-by-N matrix A.
67 * On exit, A has been destroyed.
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,M).
71 *
72 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
73 * On entry, the M-by-NRHS right hand side matrix B.
74 * On exit, B is overwritten by the N-by-NRHS solution matrix X.
75 * If m >= n and RANK = n, the residual sum-of-squares for
76 * the solution in the i-th column is given by the sum of
77 * squares of the modulus of elements n+1:m in that column.
78 *
79 * LDB (input) INTEGER
80 * The leading dimension of the array B. LDB >= max(1,M,N).
81 *
82 * S (output) DOUBLE PRECISION array, dimension (min(M,N))
83 * The singular values of A in decreasing order.
84 * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
85 *
86 * RCOND (input) DOUBLE PRECISION
87 * RCOND is used to determine the effective rank of A.
88 * Singular values S(i) <= RCOND*S(1) are treated as zero.
89 * If RCOND < 0, machine precision is used instead.
90 *
91 * RANK (output) INTEGER
92 * The effective rank of A, i.e., the number of singular values
93 * which are greater than RCOND*S(1).
94 *
95 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
96 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
97 *
98 * LWORK (input) INTEGER
99 * The dimension of the array WORK. LWORK must be at least 1.
100 * The exact minimum amount of workspace needed depends on M,
101 * N and NRHS. As long as LWORK is at least
102 * 2*N + N*NRHS
103 * if M is greater than or equal to N or
104 * 2*M + M*NRHS
105 * if M is less than N, the code will execute correctly.
106 * For good performance, LWORK should generally be larger.
107 *
108 * If LWORK = -1, then a workspace query is assumed; the routine
109 * only calculates the optimal size of the array WORK and the
110 * minimum sizes of the arrays RWORK and IWORK, and returns
111 * these values as the first entries of the WORK, RWORK and
112 * IWORK arrays, and no error message related to LWORK is issued
113 * by XERBLA.
114 *
115 * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
116 * LRWORK >=
117 * 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
118 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
119 * if M is greater than or equal to N or
120 * 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
121 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
122 * if M is less than N, the code will execute correctly.
123 * SMLSIZ is returned by ILAENV and is equal to the maximum
124 * size of the subproblems at the bottom of the computation
125 * tree (usually about 25), and
126 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
127 * On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK.
128 *
129 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
130 * LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
131 * where MINMN = MIN( M,N ).
132 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
133 *
134 * INFO (output) INTEGER
135 * = 0: successful exit
136 * < 0: if INFO = -i, the i-th argument had an illegal value.
137 * > 0: the algorithm for computing the SVD failed to converge;
138 * if INFO = i, i off-diagonal elements of an intermediate
139 * bidiagonal form did not converge to zero.
140 *
141 * Further Details
142 * ===============
143 *
144 * Based on contributions by
145 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
146 * California at Berkeley, USA
147 * Osni Marques, LBNL/NERSC, USA
148 *
149 * =====================================================================
150 *
151 * .. Parameters ..
152 DOUBLE PRECISION ZERO, ONE, TWO
153 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
154 COMPLEX*16 CZERO
155 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
156 * ..
157 * .. Local Scalars ..
158 LOGICAL LQUERY
159 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
160 $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
161 $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
162 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF,
166 $ ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR,
167 $ ZUNMLQ, ZUNMQR
168 * ..
169 * .. External Functions ..
170 INTEGER ILAENV
171 DOUBLE PRECISION DLAMCH, ZLANGE
172 EXTERNAL ILAENV, DLAMCH, ZLANGE
173 * ..
174 * .. Intrinsic Functions ..
175 INTRINSIC INT, LOG, MAX, MIN, DBLE
176 * ..
177 * .. Executable Statements ..
178 *
179 * Test the input arguments.
180 *
181 INFO = 0
182 MINMN = MIN( M, N )
183 MAXMN = MAX( M, N )
184 LQUERY = ( LWORK.EQ.-1 )
185 IF( M.LT.0 ) THEN
186 INFO = -1
187 ELSE IF( N.LT.0 ) THEN
188 INFO = -2
189 ELSE IF( NRHS.LT.0 ) THEN
190 INFO = -3
191 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
192 INFO = -5
193 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
194 INFO = -7
195 END IF
196 *
197 * Compute workspace.
198 * (Note: Comments in the code beginning "Workspace:" describe the
199 * minimal amount of workspace needed at that point in the code,
200 * as well as the preferred amount for good performance.
201 * NB refers to the optimal block size for the immediately
202 * following subroutine, as returned by ILAENV.)
203 *
204 IF( INFO.EQ.0 ) THEN
205 MINWRK = 1
206 MAXWRK = 1
207 LIWORK = 1
208 LRWORK = 1
209 IF( MINMN.GT.0 ) THEN
210 SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 )
211 MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 )
212 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) /
213 $ LOG( TWO ) ) + 1, 0 )
214 LIWORK = 3*MINMN*NLVL + 11*MINMN
215 MM = M
216 IF( M.GE.N .AND. M.GE.MNTHR ) THEN
217 *
218 * Path 1a - overdetermined, with many more rows than
219 * columns.
220 *
221 MM = N
222 MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N,
223 $ -1, -1 ) )
224 MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M,
225 $ NRHS, N, -1 ) )
226 END IF
227 IF( M.GE.N ) THEN
228 *
229 * Path 1 - overdetermined or exactly determined.
230 *
231 LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
232 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
233 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
234 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
235 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
236 $ 'QLC', MM, NRHS, N, -1 ) )
237 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
238 $ 'ZUNMBR', 'PLN', N, NRHS, N, -1 ) )
239 MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
240 MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
241 END IF
242 IF( N.GT.M ) THEN
243 LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
244 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
245 IF( N.GE.MNTHR ) THEN
246 *
247 * Path 2a - underdetermined, with many more columns
248 * than rows.
249 *
250 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
251 $ -1 )
252 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
253 $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
254 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
255 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
256 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
257 $ 'ZUNMLQ', 'LC', N, NRHS, M, -1 ) )
258 IF( NRHS.GT.1 ) THEN
259 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
260 ELSE
261 MAXWRK = MAX( MAXWRK, M*M + 2*M )
262 END IF
263 MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
264 ! XXX: Ensure the Path 2a case below is triggered. The workspace
265 ! calculation should use queries for all routines eventually.
266 MAXWRK = MAX( MAXWRK,
267 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
268 ELSE
269 *
270 * Path 2 - underdetermined.
271 *
272 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
273 $ N, -1, -1 )
274 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
275 $ 'QLC', M, NRHS, M, -1 ) )
276 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR',
277 $ 'PLN', N, NRHS, M, -1 ) )
278 MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
279 END IF
280 MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
281 END IF
282 END IF
283 MINWRK = MIN( MINWRK, MAXWRK )
284 WORK( 1 ) = MAXWRK
285 IWORK( 1 ) = LIWORK
286 RWORK( 1 ) = LRWORK
287 *
288 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
289 INFO = -12
290 END IF
291 END IF
292 *
293 IF( INFO.NE.0 ) THEN
294 CALL XERBLA( 'ZGELSD', -INFO )
295 RETURN
296 ELSE IF( LQUERY ) THEN
297 RETURN
298 END IF
299 *
300 * Quick return if possible.
301 *
302 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
303 RANK = 0
304 RETURN
305 END IF
306 *
307 * Get machine parameters.
308 *
309 EPS = DLAMCH( 'P' )
310 SFMIN = DLAMCH( 'S' )
311 SMLNUM = SFMIN / EPS
312 BIGNUM = ONE / SMLNUM
313 CALL DLABAD( SMLNUM, BIGNUM )
314 *
315 * Scale A if max entry outside range [SMLNUM,BIGNUM].
316 *
317 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
318 IASCL = 0
319 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
320 *
321 * Scale matrix norm up to SMLNUM
322 *
323 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
324 IASCL = 1
325 ELSE IF( ANRM.GT.BIGNUM ) THEN
326 *
327 * Scale matrix norm down to BIGNUM.
328 *
329 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
330 IASCL = 2
331 ELSE IF( ANRM.EQ.ZERO ) THEN
332 *
333 * Matrix all zero. Return zero solution.
334 *
335 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
336 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
337 RANK = 0
338 GO TO 10
339 END IF
340 *
341 * Scale B if max entry outside range [SMLNUM,BIGNUM].
342 *
343 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
344 IBSCL = 0
345 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
346 *
347 * Scale matrix norm up to SMLNUM.
348 *
349 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
350 IBSCL = 1
351 ELSE IF( BNRM.GT.BIGNUM ) THEN
352 *
353 * Scale matrix norm down to BIGNUM.
354 *
355 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
356 IBSCL = 2
357 END IF
358 *
359 * If M < N make sure B(M+1:N,:) = 0
360 *
361 IF( M.LT.N )
362 $ CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
363 *
364 * Overdetermined case.
365 *
366 IF( M.GE.N ) THEN
367 *
368 * Path 1 - overdetermined or exactly determined.
369 *
370 MM = M
371 IF( M.GE.MNTHR ) THEN
372 *
373 * Path 1a - overdetermined, with many more rows than columns
374 *
375 MM = N
376 ITAU = 1
377 NWORK = ITAU + N
378 *
379 * Compute A=Q*R.
380 * (RWorkspace: need N)
381 * (CWorkspace: need N, prefer N*NB)
382 *
383 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
384 $ LWORK-NWORK+1, INFO )
385 *
386 * Multiply B by transpose(Q).
387 * (RWorkspace: need N)
388 * (CWorkspace: need NRHS, prefer NRHS*NB)
389 *
390 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
391 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
392 *
393 * Zero out below R.
394 *
395 IF( N.GT.1 ) THEN
396 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
397 $ LDA )
398 END IF
399 END IF
400 *
401 ITAUQ = 1
402 ITAUP = ITAUQ + N
403 NWORK = ITAUP + N
404 IE = 1
405 NRWORK = IE + N
406 *
407 * Bidiagonalize R in A.
408 * (RWorkspace: need N)
409 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
410 *
411 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
412 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
413 $ INFO )
414 *
415 * Multiply B by transpose of left bidiagonalizing vectors of R.
416 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
417 *
418 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
419 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
420 *
421 * Solve the bidiagonal least squares problem.
422 *
423 CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
424 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
425 $ IWORK, INFO )
426 IF( INFO.NE.0 ) THEN
427 GO TO 10
428 END IF
429 *
430 * Multiply B by right bidiagonalizing vectors of R.
431 *
432 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
433 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
434 *
435 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
436 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
437 *
438 * Path 2a - underdetermined, with many more columns than rows
439 * and sufficient workspace for an efficient algorithm.
440 *
441 LDWORK = M
442 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
443 $ M*LDA+M+M*NRHS ) )LDWORK = LDA
444 ITAU = 1
445 NWORK = M + 1
446 *
447 * Compute A=L*Q.
448 * (CWorkspace: need 2*M, prefer M+M*NB)
449 *
450 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
451 $ LWORK-NWORK+1, INFO )
452 IL = NWORK
453 *
454 * Copy L to WORK(IL), zeroing out above its diagonal.
455 *
456 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
457 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
458 $ LDWORK )
459 ITAUQ = IL + LDWORK*M
460 ITAUP = ITAUQ + M
461 NWORK = ITAUP + M
462 IE = 1
463 NRWORK = IE + M
464 *
465 * Bidiagonalize L in WORK(IL).
466 * (RWorkspace: need M)
467 * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
468 *
469 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
470 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
471 $ LWORK-NWORK+1, INFO )
472 *
473 * Multiply B by transpose of left bidiagonalizing vectors of L.
474 * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
475 *
476 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
477 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
478 $ LWORK-NWORK+1, INFO )
479 *
480 * Solve the bidiagonal least squares problem.
481 *
482 CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
483 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
484 $ IWORK, INFO )
485 IF( INFO.NE.0 ) THEN
486 GO TO 10
487 END IF
488 *
489 * Multiply B by right bidiagonalizing vectors of L.
490 *
491 CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
492 $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
493 $ LWORK-NWORK+1, INFO )
494 *
495 * Zero out below first M rows of B.
496 *
497 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
498 NWORK = ITAU + M
499 *
500 * Multiply transpose(Q) by B.
501 * (CWorkspace: need NRHS, prefer NRHS*NB)
502 *
503 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
504 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
505 *
506 ELSE
507 *
508 * Path 2 - remaining underdetermined cases.
509 *
510 ITAUQ = 1
511 ITAUP = ITAUQ + M
512 NWORK = ITAUP + M
513 IE = 1
514 NRWORK = IE + M
515 *
516 * Bidiagonalize A.
517 * (RWorkspace: need M)
518 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
519 *
520 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
521 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
522 $ INFO )
523 *
524 * Multiply B by transpose of left bidiagonalizing vectors.
525 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
526 *
527 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
528 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
529 *
530 * Solve the bidiagonal least squares problem.
531 *
532 CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
533 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
534 $ IWORK, INFO )
535 IF( INFO.NE.0 ) THEN
536 GO TO 10
537 END IF
538 *
539 * Multiply B by right bidiagonalizing vectors of A.
540 *
541 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
542 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
543 *
544 END IF
545 *
546 * Undo scaling.
547 *
548 IF( IASCL.EQ.1 ) THEN
549 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
550 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
551 $ INFO )
552 ELSE IF( IASCL.EQ.2 ) THEN
553 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
554 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
555 $ INFO )
556 END IF
557 IF( IBSCL.EQ.1 ) THEN
558 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
559 ELSE IF( IBSCL.EQ.2 ) THEN
560 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
561 END IF
562 *
563 10 CONTINUE
564 WORK( 1 ) = MAXWRK
565 IWORK( 1 ) = LIWORK
566 RWORK( 1 ) = LRWORK
567 RETURN
568 *
569 * End of ZGELSD
570 *
571 END