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

CGELSS solves overdetermined or underdetermined systems for GE matrices

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

Purpose: 
 CGELSS 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 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 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 REAL 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 REAL
          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 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 REAL 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 cgelss.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  REAL rcond
189 * ..
190 * .. Array Arguments ..
191  REAL rwork( * ), s( * )
192  COMPLEX a( lda, * ), b( ldb, * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  REAL zero, one
199  parameter ( zero = 0.0e+0, one = 1.0e+0 )
200  COMPLEX czero, cone
201  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
202  $ cone = ( 1.0e+0, 0.0e+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_cgeqrf, lwork_cunmqr, lwork_cgebrd,
210  $ lwork_cunmbr, lwork_cungbr, lwork_cunmlq,
211  $ lwork_cgelqf
212  REAL anrm, bignum, bnrm, eps, sfmin, smlnum, thr
213 * ..
214 * .. Local Arrays ..
215  COMPLEX dum( 1 )
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL cbdsqr, ccopy, cgebrd, cgelqf, cgemm, cgemv,
219  $ cgeqrf, clacpy, clascl, claset, csrscl, cungbr,
220  $ cunmbr, cunmlq, cunmqr, slabad, slascl, slaset,
221  $ xerbla
222 * ..
223 * .. External Functions ..
224  INTEGER ilaenv
225  REAL clange, slamch
226  EXTERNAL ilaenv, clange, slamch
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, 'CGELSS', ' ', 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 CGEQRF
271  CALL cgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
272  lwork_cgeqrf=dum(1)
273 * Compute space needed for CUNMQR
274  CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
275  $ ldb, dum(1), -1, info )
276  lwork_cunmqr=dum(1)
277  mm = n
278  maxwrk = max( maxwrk, n + n*ilaenv( 1, 'CGEQRF', ' ', m,
279  $ n, -1, -1 ) )
280  maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'CUNMQR', '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 CGEBRD
288  CALL cgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
289  $ -1, info )
290  lwork_cgebrd=dum(1)
291 * Compute space needed for CUNMBR
292  CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
293  $ b, ldb, dum(1), -1, info )
294  lwork_cunmbr=dum(1)
295 * Compute space needed for CUNGBR
296  CALL cungbr( 'P', n, n, n, a, lda, dum(1),
297  $ dum(1), -1, info )
298  lwork_cungbr=dum(1)
299 * Compute total workspace needed
300  maxwrk = max( maxwrk, 2*n + lwork_cgebrd )
301  maxwrk = max( maxwrk, 2*n + lwork_cunmbr )
302  maxwrk = max( maxwrk, 2*n + lwork_cungbr )
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 CGELQF
314  CALL cgelqf( m, n, a, lda, dum(1), dum(1),
315  $ -1, info )
316  lwork_cgelqf=dum(1)
317 * Compute space needed for CGEBRD
318  CALL cgebrd( m, m, a, lda, s, s, dum(1), dum(1),
319  $ dum(1), -1, info )
320  lwork_cgebrd=dum(1)
321 * Compute space needed for CUNMBR
322  CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
323  $ dum(1), b, ldb, dum(1), -1, info )
324  lwork_cunmbr=dum(1)
325 * Compute space needed for CUNGBR
326  CALL cungbr( 'P', m, m, m, a, lda, dum(1),
327  $ dum(1), -1, info )
328  lwork_cungbr=dum(1)
329 * Compute space needed for CUNMLQ
330  CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
331  $ b, ldb, dum(1), -1, info )
332  lwork_cunmlq=dum(1)
333 * Compute total workspace needed
334  maxwrk = m + lwork_cgelqf
335  maxwrk = max( maxwrk, 3*m + m*m + lwork_cgebrd )
336  maxwrk = max( maxwrk, 3*m + m*m + lwork_cunmbr )
337  maxwrk = max( maxwrk, 3*m + m*m + lwork_cungbr )
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_cunmlq )
344  ELSE
345 *
346 * Path 2 - underdetermined
347 *
348 * Compute space needed for CGEBRD
349  CALL cgebrd( m, n, a, lda, s, s, dum(1), dum(1),
350  $ dum(1), -1, info )
351  lwork_cgebrd=dum(1)
352 * Compute space needed for CUNMBR
353  CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
354  $ dum(1), b, ldb, dum(1), -1, info )
355  lwork_cunmbr=dum(1)
356 * Compute space needed for CUNGBR
357  CALL cungbr( 'P', m, n, m, a, lda, dum(1),
358  $ dum(1), -1, info )
359  lwork_cungbr=dum(1)
360  maxwrk = 2*m + lwork_cgebrd
361  maxwrk = max( maxwrk, 2*m + lwork_cunmbr )
362  maxwrk = max( maxwrk, 2*m + lwork_cungbr )
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( 'CGELSS', -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 = slamch( 'P' )
391  sfmin = slamch( 'S' )
392  smlnum = sfmin / eps
393  bignum = one / smlnum
394  CALL slabad( smlnum, bignum )
395 *
396 * Scale A if max element outside range [SMLNUM,BIGNUM]
397 *
398  anrm = clange( '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 clascl( '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 clascl( '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 claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
417  CALL slaset( '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 = clange( '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 clascl( '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 clascl( '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 cgeqrf( 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 cunmqr( '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 claset( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
524  rank = rank + 1
525  ELSE
526  CALL claset( '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 cgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
536  $ czero, work, ldb )
537  CALL clacpy( '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 cgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
543  $ ldb, czero, work, n )
544  CALL clacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
545  20 CONTINUE
546  ELSE
547  CALL cgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
548  CALL ccopy( 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 cgelqf( 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 clacpy( 'L', m, m, a, lda, work( il ), ldwork )
576  CALL claset( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
627  rank = rank + 1
628  ELSE
629  CALL claset( '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 cgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
640  $ b, ldb, czero, work( iwork ), ldb )
641  CALL clacpy( '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 cgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
647  $ b( 1, i ), ldb, czero, work( iwork ), m )
648  CALL clacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
649  $ ldb )
650  40 CONTINUE
651  ELSE
652  CALL cgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
653  $ 1, czero, work( iwork ), 1 )
654  CALL ccopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
655  END IF
656 *
657 * Zero out below first M rows of B
658 *
659  CALL claset( '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 cunmlq( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
721  rank = rank + 1
722  ELSE
723  CALL claset( '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 cgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
733  $ czero, work, ldb )
734  CALL clacpy( '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 cgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
740  $ ldb, czero, work, n )
741  CALL clacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
742  60 CONTINUE
743  ELSE
744  CALL cgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
745  CALL ccopy( n, work, 1, b, 1 )
746  END IF
747  END IF
748 *
749 * Undo scaling
750 *
751  IF( iascl.EQ.1 ) THEN
752  CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
753  CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
754  $ info )
755  ELSE IF( iascl.EQ.2 ) THEN
756  CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
757  CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
758  $ info )
759  END IF
760  IF( ibscl.EQ.1 ) THEN
761  CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
762  ELSE IF( ibscl.EQ.2 ) THEN
763  CALL clascl( '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 CGELSS
770 *
cungbr
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
Definition: cungbr.f:159
cunmbr
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:199
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
cbdsqr
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
Definition: cbdsqr.f:224
clange
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cunmqr
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
cgelqf
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:137
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:138
cgebrd
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
Definition: cgebrd.f:208
csrscl
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: csrscl.f:86
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
cunmlq
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:170
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189

Here is the call graph for this function:

Here is the caller graph for this function: