LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zgelss ( integer  M,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  S,
double precision  RCOND,
integer  RANK,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGELSS solves overdetermined or underdetermined systems for GE matrices

Download ZGELSS + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 ZGELSS computes the minimum norm solution to a complex linear
 least squares problem:

 Minimize 2-norm(| b - A*x |).

 using the singular value decomposition (SVD) of A. A is an M-by-N
 matrix which may be rank-deficient.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
 X.

 The effective rank of A is determined by treating as zero those
 singular values which are less than RCOND times the largest singular
 value.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A. N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X. NRHS >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the first min(m,n) rows of A are overwritten with
          its right singular vectors, stored rowwise.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          On entry, the M-by-NRHS right hand side matrix B.
          On exit, B is overwritten by the N-by-NRHS solution matrix X.
          If m >= n and RANK = n, the residual sum-of-squares for
          the solution in the i-th column is given by the sum of
          squares of the modulus of elements n+1:m in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,M,N).
[out]S
          S is DOUBLE PRECISION array, dimension (min(M,N))
          The singular values of A in decreasing order.
          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
[in]RCOND
          RCOND is DOUBLE PRECISION
          RCOND is used to determine the effective rank of A.
          Singular values S(i) <= RCOND*S(1) are treated as zero.
          If RCOND < 0, machine precision is used instead.
[out]RANK
          RANK is INTEGER
          The effective rank of A, i.e., the number of singular values
          which are greater than RCOND*S(1).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= 1, and also:
          LWORK >=  2*min(M,N) + max(M,N,NRHS)
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (5*min(M,N))
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  the algorithm for computing the SVD failed to converge;
                if INFO = i, i off-diagonal elements of an intermediate
                bidiagonal form did not converge to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 180 of file zgelss.f.

180 *
181 * -- LAPACK driver routine (version 3.6.1) --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 * June 2016
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,
219  $ zdrscl, zgebrd, zgelqf, zgemm, zgemv, zgeqrf,
220  $ zlacpy, zlascl, zlaset, zungbr, zunmbr, zunmlq,
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, s, dum(1), dum(1), dum(1),
289  $ -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, s, dum(1), dum(1),
319  $ 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, s, dum(1), dum(1),
350  $ 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 *
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
zdrscl
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zungbr
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:159
zgeqrf
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zunmbr
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:198
dlascl
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
zgebrd
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:207
zunmqr
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
zlange
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
zunmlq
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:169
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
zbdsqr
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
Definition: zbdsqr.f:224
zgelqf
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:137
zlascl
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145

Here is the call graph for this function:

Here is the caller graph for this function: