LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zgelss()

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.

Definition at line 176 of file zgelss.f.

178 *
179 * -- LAPACK driver routine --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 *
183 * .. Scalar Arguments ..
184  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
185  DOUBLE PRECISION RCOND
186 * ..
187 * .. Array Arguments ..
188  DOUBLE PRECISION RWORK( * ), S( * )
189  COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  DOUBLE PRECISION ZERO, ONE
196  parameter( zero = 0.0d+0, one = 1.0d+0 )
197  COMPLEX*16 CZERO, CONE
198  parameter( czero = ( 0.0d+0, 0.0d+0 ),
199  $ cone = ( 1.0d+0, 0.0d+0 ) )
200 * ..
201 * .. Local Scalars ..
202  LOGICAL LQUERY
203  INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
204  $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
205  $ MAXWRK, MINMN, MINWRK, MM, MNTHR
206  INTEGER LWORK_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD,
207  $ LWORK_ZUNMBR, LWORK_ZUNGBR, LWORK_ZUNMLQ,
208  $ LWORK_ZGELQF
209  DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
210 * ..
211 * .. Local Arrays ..
212  COMPLEX*16 DUM( 1 )
213 * ..
214 * .. External Subroutines ..
215  EXTERNAL dlabad, dlascl, dlaset, xerbla, zbdsqr, zcopy,
218  $ zunmqr
219 * ..
220 * .. External Functions ..
221  INTEGER ILAENV
222  DOUBLE PRECISION DLAMCH, ZLANGE
223  EXTERNAL ilaenv, dlamch, zlange
224 * ..
225 * .. Intrinsic Functions ..
226  INTRINSIC max, min
227 * ..
228 * .. Executable Statements ..
229 *
230 * Test the input arguments
231 *
232  info = 0
233  minmn = min( m, n )
234  maxmn = max( m, n )
235  lquery = ( lwork.EQ.-1 )
236  IF( m.LT.0 ) THEN
237  info = -1
238  ELSE IF( n.LT.0 ) THEN
239  info = -2
240  ELSE IF( nrhs.LT.0 ) THEN
241  info = -3
242  ELSE IF( lda.LT.max( 1, m ) ) THEN
243  info = -5
244  ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
245  info = -7
246  END IF
247 *
248 * Compute workspace
249 * (Note: Comments in the code beginning "Workspace:" describe the
250 * minimal amount of workspace needed at that point in the code,
251 * as well as the preferred amount for good performance.
252 * CWorkspace refers to complex workspace, and RWorkspace refers
253 * to real workspace. NB refers to the optimal block size for the
254 * immediately following subroutine, as returned by ILAENV.)
255 *
256  IF( info.EQ.0 ) THEN
257  minwrk = 1
258  maxwrk = 1
259  IF( minmn.GT.0 ) THEN
260  mm = m
261  mnthr = ilaenv( 6, 'ZGELSS', ' ', m, n, nrhs, -1 )
262  IF( m.GE.n .AND. m.GE.mnthr ) THEN
263 *
264 * Path 1a - overdetermined, with many more rows than
265 * columns
266 *
267 * Compute space needed for ZGEQRF
268  CALL zgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
269  lwork_zgeqrf = dble( dum(1) )
270 * Compute space needed for ZUNMQR
271  CALL zunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
272  $ ldb, dum(1), -1, info )
273  lwork_zunmqr = dble( dum(1) )
274  mm = n
275  maxwrk = max( maxwrk, n + n*ilaenv( 1, 'ZGEQRF', ' ', m,
276  $ n, -1, -1 ) )
277  maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'ZUNMQR', 'LC',
278  $ m, nrhs, n, -1 ) )
279  END IF
280  IF( m.GE.n ) THEN
281 *
282 * Path 1 - overdetermined or exactly determined
283 *
284 * Compute space needed for ZGEBRD
285  CALL zgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
286  $ -1, info )
287  lwork_zgebrd = dble( dum(1) )
288 * Compute space needed for ZUNMBR
289  CALL zunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
290  $ b, ldb, dum(1), -1, info )
291  lwork_zunmbr = dble( dum(1) )
292 * Compute space needed for ZUNGBR
293  CALL zungbr( 'P', n, n, n, a, lda, dum(1),
294  $ dum(1), -1, info )
295  lwork_zungbr = dble( dum(1) )
296 * Compute total workspace needed
297  maxwrk = max( maxwrk, 2*n + lwork_zgebrd )
298  maxwrk = max( maxwrk, 2*n + lwork_zunmbr )
299  maxwrk = max( maxwrk, 2*n + lwork_zungbr )
300  maxwrk = max( maxwrk, n*nrhs )
301  minwrk = 2*n + max( nrhs, m )
302  END IF
303  IF( n.GT.m ) THEN
304  minwrk = 2*m + max( nrhs, n )
305  IF( n.GE.mnthr ) THEN
306 *
307 * Path 2a - underdetermined, with many more columns
308 * than rows
309 *
310 * Compute space needed for ZGELQF
311  CALL zgelqf( m, n, a, lda, dum(1), dum(1),
312  $ -1, info )
313  lwork_zgelqf = dble( dum(1) )
314 * Compute space needed for ZGEBRD
315  CALL zgebrd( m, m, a, lda, s, s, dum(1), dum(1),
316  $ dum(1), -1, info )
317  lwork_zgebrd = dble( dum(1) )
318 * Compute space needed for ZUNMBR
319  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
320  $ dum(1), b, ldb, dum(1), -1, info )
321  lwork_zunmbr = dble( dum(1) )
322 * Compute space needed for ZUNGBR
323  CALL zungbr( 'P', m, m, m, a, lda, dum(1),
324  $ dum(1), -1, info )
325  lwork_zungbr = dble( dum(1) )
326 * Compute space needed for ZUNMLQ
327  CALL zunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
328  $ b, ldb, dum(1), -1, info )
329  lwork_zunmlq = dble( dum(1) )
330 * Compute total workspace needed
331  maxwrk = m + lwork_zgelqf
332  maxwrk = max( maxwrk, 3*m + m*m + lwork_zgebrd )
333  maxwrk = max( maxwrk, 3*m + m*m + lwork_zunmbr )
334  maxwrk = max( maxwrk, 3*m + m*m + lwork_zungbr )
335  IF( nrhs.GT.1 ) THEN
336  maxwrk = max( maxwrk, m*m + m + m*nrhs )
337  ELSE
338  maxwrk = max( maxwrk, m*m + 2*m )
339  END IF
340  maxwrk = max( maxwrk, m + lwork_zunmlq )
341  ELSE
342 *
343 * Path 2 - underdetermined
344 *
345 * Compute space needed for ZGEBRD
346  CALL zgebrd( m, n, a, lda, s, s, dum(1), dum(1),
347  $ dum(1), -1, info )
348  lwork_zgebrd = dble( dum(1) )
349 * Compute space needed for ZUNMBR
350  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
351  $ dum(1), b, ldb, dum(1), -1, info )
352  lwork_zunmbr = dble( dum(1) )
353 * Compute space needed for ZUNGBR
354  CALL zungbr( 'P', m, n, m, a, lda, dum(1),
355  $ dum(1), -1, info )
356  lwork_zungbr = dble( dum(1) )
357  maxwrk = 2*m + lwork_zgebrd
358  maxwrk = max( maxwrk, 2*m + lwork_zunmbr )
359  maxwrk = max( maxwrk, 2*m + lwork_zungbr )
360  maxwrk = max( maxwrk, n*nrhs )
361  END IF
362  END IF
363  maxwrk = max( minwrk, maxwrk )
364  END IF
365  work( 1 ) = maxwrk
366 *
367  IF( lwork.LT.minwrk .AND. .NOT.lquery )
368  $ info = -12
369  END IF
370 *
371  IF( info.NE.0 ) THEN
372  CALL xerbla( 'ZGELSS', -info )
373  RETURN
374  ELSE IF( lquery ) THEN
375  RETURN
376  END IF
377 *
378 * Quick return if possible
379 *
380  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
381  rank = 0
382  RETURN
383  END IF
384 *
385 * Get machine parameters
386 *
387  eps = dlamch( 'P' )
388  sfmin = dlamch( 'S' )
389  smlnum = sfmin / eps
390  bignum = one / smlnum
391  CALL dlabad( smlnum, bignum )
392 *
393 * Scale A if max element outside range [SMLNUM,BIGNUM]
394 *
395  anrm = zlange( 'M', m, n, a, lda, rwork )
396  iascl = 0
397  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
398 *
399 * Scale matrix norm up to SMLNUM
400 *
401  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
402  iascl = 1
403  ELSE IF( anrm.GT.bignum ) THEN
404 *
405 * Scale matrix norm down to BIGNUM
406 *
407  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
408  iascl = 2
409  ELSE IF( anrm.EQ.zero ) THEN
410 *
411 * Matrix all zero. Return zero solution.
412 *
413  CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
414  CALL dlaset( 'F', minmn, 1, zero, zero, s, minmn )
415  rank = 0
416  GO TO 70
417  END IF
418 *
419 * Scale B if max element outside range [SMLNUM,BIGNUM]
420 *
421  bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
422  ibscl = 0
423  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
424 *
425 * Scale matrix norm up to SMLNUM
426 *
427  CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
428  ibscl = 1
429  ELSE IF( bnrm.GT.bignum ) THEN
430 *
431 * Scale matrix norm down to BIGNUM
432 *
433  CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
434  ibscl = 2
435  END IF
436 *
437 * Overdetermined case
438 *
439  IF( m.GE.n ) THEN
440 *
441 * Path 1 - overdetermined or exactly determined
442 *
443  mm = m
444  IF( m.GE.mnthr ) THEN
445 *
446 * Path 1a - overdetermined, with many more rows than columns
447 *
448  mm = n
449  itau = 1
450  iwork = itau + n
451 *
452 * Compute A=Q*R
453 * (CWorkspace: need 2*N, prefer N+N*NB)
454 * (RWorkspace: none)
455 *
456  CALL zgeqrf( m, n, a, lda, work( itau ), work( iwork ),
457  $ lwork-iwork+1, info )
458 *
459 * Multiply B by transpose(Q)
460 * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
461 * (RWorkspace: none)
462 *
463  CALL zunmqr( 'L', 'C', m, nrhs, n, a, lda, work( itau ), b,
464  $ ldb, work( iwork ), lwork-iwork+1, info )
465 *
466 * Zero out below R
467 *
468  IF( n.GT.1 )
469  $ CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
470  $ lda )
471  END IF
472 *
473  ie = 1
474  itauq = 1
475  itaup = itauq + n
476  iwork = itaup + n
477 *
478 * Bidiagonalize R in A
479 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
480 * (RWorkspace: need N)
481 *
482  CALL zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
483  $ work( itaup ), work( iwork ), lwork-iwork+1,
484  $ info )
485 *
486 * Multiply B by transpose of left bidiagonalizing vectors of R
487 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
488 * (RWorkspace: none)
489 *
490  CALL zunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, work( itauq ),
491  $ b, ldb, work( iwork ), lwork-iwork+1, info )
492 *
493 * Generate right bidiagonalizing vectors of R in A
494 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
495 * (RWorkspace: none)
496 *
497  CALL zungbr( 'P', n, n, n, a, lda, work( itaup ),
498  $ work( iwork ), lwork-iwork+1, info )
499  irwork = ie + n
500 *
501 * Perform bidiagonal QR iteration
502 * multiply B by transpose of left singular vectors
503 * compute right singular vectors in A
504 * (CWorkspace: none)
505 * (RWorkspace: need BDSPAC)
506 *
507  CALL zbdsqr( 'U', n, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
508  $ 1, b, ldb, rwork( irwork ), info )
509  IF( info.NE.0 )
510  $ GO TO 70
511 *
512 * Multiply B by reciprocals of singular values
513 *
514  thr = max( rcond*s( 1 ), sfmin )
515  IF( rcond.LT.zero )
516  $ thr = max( eps*s( 1 ), sfmin )
517  rank = 0
518  DO 10 i = 1, n
519  IF( s( i ).GT.thr ) THEN
520  CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
521  rank = rank + 1
522  ELSE
523  CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
524  END IF
525  10 CONTINUE
526 *
527 * Multiply B by right singular vectors
528 * (CWorkspace: need N, prefer N*NRHS)
529 * (RWorkspace: none)
530 *
531  IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
532  CALL zgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
533  $ czero, work, ldb )
534  CALL zlacpy( 'G', n, nrhs, work, ldb, b, ldb )
535  ELSE IF( nrhs.GT.1 ) THEN
536  chunk = lwork / n
537  DO 20 i = 1, nrhs, chunk
538  bl = min( nrhs-i+1, chunk )
539  CALL zgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
540  $ ldb, czero, work, n )
541  CALL zlacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
542  20 CONTINUE
543  ELSE
544  CALL zgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
545  CALL zcopy( n, work, 1, b, 1 )
546  END IF
547 *
548  ELSE IF( n.GE.mnthr .AND. lwork.GE.3*m+m*m+max( m, nrhs, n-2*m ) )
549  $ THEN
550 *
551 * Underdetermined case, M much less than N
552 *
553 * Path 2a - underdetermined, with many more columns than rows
554 * and sufficient workspace for an efficient algorithm
555 *
556  ldwork = m
557  IF( lwork.GE.3*m+m*lda+max( m, nrhs, n-2*m ) )
558  $ ldwork = lda
559  itau = 1
560  iwork = m + 1
561 *
562 * Compute A=L*Q
563 * (CWorkspace: need 2*M, prefer M+M*NB)
564 * (RWorkspace: none)
565 *
566  CALL zgelqf( m, n, a, lda, work( itau ), work( iwork ),
567  $ lwork-iwork+1, info )
568  il = iwork
569 *
570 * Copy L to WORK(IL), zeroing out above it
571 *
572  CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwork )
573  CALL zlaset( 'U', m-1, m-1, czero, czero, work( il+ldwork ),
574  $ ldwork )
575  ie = 1
576  itauq = il + ldwork*m
577  itaup = itauq + m
578  iwork = itaup + m
579 *
580 * Bidiagonalize L in WORK(IL)
581 * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
582 * (RWorkspace: need M)
583 *
584  CALL zgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
585  $ work( itauq ), work( itaup ), work( iwork ),
586  $ lwork-iwork+1, info )
587 *
588 * Multiply B by transpose of left bidiagonalizing vectors of L
589 * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
590 * (RWorkspace: none)
591 *
592  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, m, work( il ), ldwork,
593  $ work( itauq ), b, ldb, work( iwork ),
594  $ lwork-iwork+1, info )
595 *
596 * Generate right bidiagonalizing vectors of R in WORK(IL)
597 * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
598 * (RWorkspace: none)
599 *
600  CALL zungbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
601  $ work( iwork ), lwork-iwork+1, info )
602  irwork = ie + m
603 *
604 * Perform bidiagonal QR iteration, computing right singular
605 * vectors of L in WORK(IL) and multiplying B by transpose of
606 * left singular vectors
607 * (CWorkspace: need M*M)
608 * (RWorkspace: need BDSPAC)
609 *
610  CALL zbdsqr( 'U', m, m, 0, nrhs, s, rwork( ie ), work( il ),
611  $ ldwork, a, lda, b, ldb, rwork( irwork ), info )
612  IF( info.NE.0 )
613  $ GO TO 70
614 *
615 * Multiply B by reciprocals of singular values
616 *
617  thr = max( rcond*s( 1 ), sfmin )
618  IF( rcond.LT.zero )
619  $ thr = max( eps*s( 1 ), sfmin )
620  rank = 0
621  DO 30 i = 1, m
622  IF( s( i ).GT.thr ) THEN
623  CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
624  rank = rank + 1
625  ELSE
626  CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
627  END IF
628  30 CONTINUE
629  iwork = il + m*ldwork
630 *
631 * Multiply B by right singular vectors of L in WORK(IL)
632 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
633 * (RWorkspace: none)
634 *
635  IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
636  CALL zgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
637  $ b, ldb, czero, work( iwork ), ldb )
638  CALL zlacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
639  ELSE IF( nrhs.GT.1 ) THEN
640  chunk = ( lwork-iwork+1 ) / m
641  DO 40 i = 1, nrhs, chunk
642  bl = min( nrhs-i+1, chunk )
643  CALL zgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
644  $ b( 1, i ), ldb, czero, work( iwork ), m )
645  CALL zlacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
646  $ ldb )
647  40 CONTINUE
648  ELSE
649  CALL zgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
650  $ 1, czero, work( iwork ), 1 )
651  CALL zcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
652  END IF
653 *
654 * Zero out below first M rows of B
655 *
656  CALL zlaset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
657  iwork = itau + m
658 *
659 * Multiply transpose(Q) by B
660 * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
661 * (RWorkspace: none)
662 *
663  CALL zunmlq( 'L', 'C', n, nrhs, m, a, lda, work( itau ), b,
664  $ ldb, work( iwork ), lwork-iwork+1, info )
665 *
666  ELSE
667 *
668 * Path 2 - remaining underdetermined cases
669 *
670  ie = 1
671  itauq = 1
672  itaup = itauq + m
673  iwork = itaup + m
674 *
675 * Bidiagonalize A
676 * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
677 * (RWorkspace: need N)
678 *
679  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
680  $ work( itaup ), work( iwork ), lwork-iwork+1,
681  $ info )
682 *
683 * Multiply B by transpose of left bidiagonalizing vectors
684 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
685 * (RWorkspace: none)
686 *
687  CALL zunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda, work( itauq ),
688  $ b, ldb, work( iwork ), lwork-iwork+1, info )
689 *
690 * Generate right bidiagonalizing vectors in A
691 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
692 * (RWorkspace: none)
693 *
694  CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
695  $ work( iwork ), lwork-iwork+1, info )
696  irwork = ie + m
697 *
698 * Perform bidiagonal QR iteration,
699 * computing right singular vectors of A in A and
700 * multiplying B by transpose of left singular vectors
701 * (CWorkspace: none)
702 * (RWorkspace: need BDSPAC)
703 *
704  CALL zbdsqr( 'L', m, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
705  $ 1, b, ldb, rwork( irwork ), info )
706  IF( info.NE.0 )
707  $ GO TO 70
708 *
709 * Multiply B by reciprocals of singular values
710 *
711  thr = max( rcond*s( 1 ), sfmin )
712  IF( rcond.LT.zero )
713  $ thr = max( eps*s( 1 ), sfmin )
714  rank = 0
715  DO 50 i = 1, m
716  IF( s( i ).GT.thr ) THEN
717  CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
718  rank = rank + 1
719  ELSE
720  CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
721  END IF
722  50 CONTINUE
723 *
724 * Multiply B by right singular vectors of A
725 * (CWorkspace: need N, prefer N*NRHS)
726 * (RWorkspace: none)
727 *
728  IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
729  CALL zgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
730  $ czero, work, ldb )
731  CALL zlacpy( 'G', n, nrhs, work, ldb, b, ldb )
732  ELSE IF( nrhs.GT.1 ) THEN
733  chunk = lwork / n
734  DO 60 i = 1, nrhs, chunk
735  bl = min( nrhs-i+1, chunk )
736  CALL zgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
737  $ ldb, czero, work, n )
738  CALL zlacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
739  60 CONTINUE
740  ELSE
741  CALL zgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
742  CALL zcopy( n, work, 1, b, 1 )
743  END IF
744  END IF
745 *
746 * Undo scaling
747 *
748  IF( iascl.EQ.1 ) THEN
749  CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
750  CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
751  $ info )
752  ELSE IF( iascl.EQ.2 ) THEN
753  CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
754  CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
755  $ info )
756  END IF
757  IF( ibscl.EQ.1 ) THEN
758  CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
759  ELSE IF( ibscl.EQ.2 ) THEN
760  CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
761  END IF
762  70 CONTINUE
763  work( 1 ) = maxwrk
764  RETURN
765 *
766 * End of ZGELSS
767 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:143
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:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:157
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:115
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:143
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:205
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:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
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:106
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:84
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:167
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
Definition: zbdsqr.f:222
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:196
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
Here is the call graph for this function:
Here is the caller graph for this function: