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