LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cgelss()

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.

Definition at line 176 of file cgelss.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 REAL RCOND
186* ..
187* .. Array Arguments ..
188 REAL RWORK( * ), S( * )
189 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ZERO, ONE
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197 COMPLEX CZERO, CONE
198 parameter( czero = ( 0.0e+0, 0.0e+0 ),
199 $ cone = ( 1.0e+0, 0.0e+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_CGEQRF, LWORK_CUNMQR, LWORK_CGEBRD,
207 $ LWORK_CUNMBR, LWORK_CUNGBR, LWORK_CUNMLQ,
208 $ LWORK_CGELQF
209 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
210* ..
211* .. Local Arrays ..
212 COMPLEX DUM( 1 )
213* ..
214* .. External Subroutines ..
215 EXTERNAL cbdsqr, ccopy, cgebrd, cgelqf, cgemm, cgemv,
216 $ cgeqrf, clacpy, clascl, claset, csrscl, cungbr,
217 $ cunmbr, cunmlq, cunmqr, slascl, slaset, xerbla
218* ..
219* .. External Functions ..
220 INTEGER ILAENV
221 REAL CLANGE, SLAMCH, SROUNDUP_LWORK
222 EXTERNAL ilaenv, clange, slamch, sroundup_lwork
223* ..
224* .. Intrinsic Functions ..
225 INTRINSIC max, min
226* ..
227* .. Executable Statements ..
228*
229* Test the input arguments
230*
231 info = 0
232 minmn = min( m, n )
233 maxmn = max( m, n )
234 lquery = ( lwork.EQ.-1 )
235 IF( m.LT.0 ) THEN
236 info = -1
237 ELSE IF( n.LT.0 ) THEN
238 info = -2
239 ELSE IF( nrhs.LT.0 ) THEN
240 info = -3
241 ELSE IF( lda.LT.max( 1, m ) ) THEN
242 info = -5
243 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
244 info = -7
245 END IF
246*
247* Compute workspace
248* (Note: Comments in the code beginning "Workspace:" describe the
249* minimal amount of workspace needed at that point in the code,
250* as well as the preferred amount for good performance.
251* CWorkspace refers to complex workspace, and RWorkspace refers
252* to real workspace. NB refers to the optimal block size for the
253* immediately following subroutine, as returned by ILAENV.)
254*
255 IF( info.EQ.0 ) THEN
256 minwrk = 1
257 maxwrk = 1
258 IF( minmn.GT.0 ) THEN
259 mm = m
260 mnthr = ilaenv( 6, 'CGELSS', ' ', m, n, nrhs, -1 )
261 IF( m.GE.n .AND. m.GE.mnthr ) THEN
262*
263* Path 1a - overdetermined, with many more rows than
264* columns
265*
266* Compute space needed for CGEQRF
267 CALL cgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
268 lwork_cgeqrf = int( dum(1) )
269* Compute space needed for CUNMQR
270 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
271 $ ldb, dum(1), -1, info )
272 lwork_cunmqr = int( dum(1) )
273 mm = n
274 maxwrk = max( maxwrk, n + n*ilaenv( 1, 'CGEQRF', ' ', m,
275 $ n, -1, -1 ) )
276 maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'CUNMQR', 'LC',
277 $ m, nrhs, n, -1 ) )
278 END IF
279 IF( m.GE.n ) THEN
280*
281* Path 1 - overdetermined or exactly determined
282*
283* Compute space needed for CGEBRD
284 CALL cgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
285 $ -1, info )
286 lwork_cgebrd = int( dum(1) )
287* Compute space needed for CUNMBR
288 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
289 $ b, ldb, dum(1), -1, info )
290 lwork_cunmbr = int( dum(1) )
291* Compute space needed for CUNGBR
292 CALL cungbr( 'P', n, n, n, a, lda, dum(1),
293 $ dum(1), -1, info )
294 lwork_cungbr = int( dum(1) )
295* Compute total workspace needed
296 maxwrk = max( maxwrk, 2*n + lwork_cgebrd )
297 maxwrk = max( maxwrk, 2*n + lwork_cunmbr )
298 maxwrk = max( maxwrk, 2*n + lwork_cungbr )
299 maxwrk = max( maxwrk, n*nrhs )
300 minwrk = 2*n + max( nrhs, m )
301 END IF
302 IF( n.GT.m ) THEN
303 minwrk = 2*m + max( nrhs, n )
304 IF( n.GE.mnthr ) THEN
305*
306* Path 2a - underdetermined, with many more columns
307* than rows
308*
309* Compute space needed for CGELQF
310 CALL cgelqf( m, n, a, lda, dum(1), dum(1),
311 $ -1, info )
312 lwork_cgelqf = int( dum(1) )
313* Compute space needed for CGEBRD
314 CALL cgebrd( m, m, a, lda, s, s, dum(1), dum(1),
315 $ dum(1), -1, info )
316 lwork_cgebrd = int( dum(1) )
317* Compute space needed for CUNMBR
318 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
319 $ dum(1), b, ldb, dum(1), -1, info )
320 lwork_cunmbr = int( dum(1) )
321* Compute space needed for CUNGBR
322 CALL cungbr( 'P', m, m, m, a, lda, dum(1),
323 $ dum(1), -1, info )
324 lwork_cungbr = int( dum(1) )
325* Compute space needed for CUNMLQ
326 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
327 $ b, ldb, dum(1), -1, info )
328 lwork_cunmlq = int( dum(1) )
329* Compute total workspace needed
330 maxwrk = m + lwork_cgelqf
331 maxwrk = max( maxwrk, 3*m + m*m + lwork_cgebrd )
332 maxwrk = max( maxwrk, 3*m + m*m + lwork_cunmbr )
333 maxwrk = max( maxwrk, 3*m + m*m + lwork_cungbr )
334 IF( nrhs.GT.1 ) THEN
335 maxwrk = max( maxwrk, m*m + m + m*nrhs )
336 ELSE
337 maxwrk = max( maxwrk, m*m + 2*m )
338 END IF
339 maxwrk = max( maxwrk, m + lwork_cunmlq )
340 ELSE
341*
342* Path 2 - underdetermined
343*
344* Compute space needed for CGEBRD
345 CALL cgebrd( m, n, a, lda, s, s, dum(1), dum(1),
346 $ dum(1), -1, info )
347 lwork_cgebrd = int( dum(1) )
348* Compute space needed for CUNMBR
349 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
350 $ dum(1), b, ldb, dum(1), -1, info )
351 lwork_cunmbr = int( dum(1) )
352* Compute space needed for CUNGBR
353 CALL cungbr( 'P', m, n, m, a, lda, dum(1),
354 $ dum(1), -1, info )
355 lwork_cungbr = int( dum(1) )
356 maxwrk = 2*m + lwork_cgebrd
357 maxwrk = max( maxwrk, 2*m + lwork_cunmbr )
358 maxwrk = max( maxwrk, 2*m + lwork_cungbr )
359 maxwrk = max( maxwrk, n*nrhs )
360 END IF
361 END IF
362 maxwrk = max( minwrk, maxwrk )
363 END IF
364 work( 1 ) = sroundup_lwork(maxwrk)
365*
366 IF( lwork.LT.minwrk .AND. .NOT.lquery )
367 $ info = -12
368 END IF
369*
370 IF( info.NE.0 ) THEN
371 CALL xerbla( 'CGELSS', -info )
372 RETURN
373 ELSE IF( lquery ) THEN
374 RETURN
375 END IF
376*
377* Quick return if possible
378*
379 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
380 rank = 0
381 RETURN
382 END IF
383*
384* Get machine parameters
385*
386 eps = slamch( 'P' )
387 sfmin = slamch( 'S' )
388 smlnum = sfmin / eps
389 bignum = one / smlnum
390*
391* Scale A if max element outside range [SMLNUM,BIGNUM]
392*
393 anrm = clange( 'M', m, n, a, lda, rwork )
394 iascl = 0
395 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
396*
397* Scale matrix norm up to SMLNUM
398*
399 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
400 iascl = 1
401 ELSE IF( anrm.GT.bignum ) THEN
402*
403* Scale matrix norm down to BIGNUM
404*
405 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
406 iascl = 2
407 ELSE IF( anrm.EQ.zero ) THEN
408*
409* Matrix all zero. Return zero solution.
410*
411 CALL claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
412 CALL slaset( 'F', minmn, 1, zero, zero, s, minmn )
413 rank = 0
414 GO TO 70
415 END IF
416*
417* Scale B if max element outside range [SMLNUM,BIGNUM]
418*
419 bnrm = clange( 'M', m, nrhs, b, ldb, rwork )
420 ibscl = 0
421 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
422*
423* Scale matrix norm up to SMLNUM
424*
425 CALL clascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
426 ibscl = 1
427 ELSE IF( bnrm.GT.bignum ) THEN
428*
429* Scale matrix norm down to BIGNUM
430*
431 CALL clascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
432 ibscl = 2
433 END IF
434*
435* Overdetermined case
436*
437 IF( m.GE.n ) THEN
438*
439* Path 1 - overdetermined or exactly determined
440*
441 mm = m
442 IF( m.GE.mnthr ) THEN
443*
444* Path 1a - overdetermined, with many more rows than columns
445*
446 mm = n
447 itau = 1
448 iwork = itau + n
449*
450* Compute A=Q*R
451* (CWorkspace: need 2*N, prefer N+N*NB)
452* (RWorkspace: none)
453*
454 CALL cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
455 $ lwork-iwork+1, info )
456*
457* Multiply B by transpose(Q)
458* (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
459* (RWorkspace: none)
460*
461 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, work( itau ), b,
462 $ ldb, work( iwork ), lwork-iwork+1, info )
463*
464* Zero out below R
465*
466 IF( n.GT.1 )
467 $ CALL claset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
468 $ lda )
469 END IF
470*
471 ie = 1
472 itauq = 1
473 itaup = itauq + n
474 iwork = itaup + n
475*
476* Bidiagonalize R in A
477* (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
478* (RWorkspace: need N)
479*
480 CALL cgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
481 $ work( itaup ), work( iwork ), lwork-iwork+1,
482 $ info )
483*
484* Multiply B by transpose of left bidiagonalizing vectors of R
485* (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
486* (RWorkspace: none)
487*
488 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, work( itauq ),
489 $ b, ldb, work( iwork ), lwork-iwork+1, info )
490*
491* Generate right bidiagonalizing vectors of R in A
492* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
493* (RWorkspace: none)
494*
495 CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
496 $ work( iwork ), lwork-iwork+1, info )
497 irwork = ie + n
498*
499* Perform bidiagonal QR iteration
500* multiply B by transpose of left singular vectors
501* compute right singular vectors in A
502* (CWorkspace: none)
503* (RWorkspace: need BDSPAC)
504*
505 CALL cbdsqr( 'U', n, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
506 $ 1, b, ldb, rwork( irwork ), info )
507 IF( info.NE.0 )
508 $ GO TO 70
509*
510* Multiply B by reciprocals of singular values
511*
512 thr = max( rcond*s( 1 ), sfmin )
513 IF( rcond.LT.zero )
514 $ thr = max( eps*s( 1 ), sfmin )
515 rank = 0
516 DO 10 i = 1, n
517 IF( s( i ).GT.thr ) THEN
518 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
519 rank = rank + 1
520 ELSE
521 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
522 END IF
523 10 CONTINUE
524*
525* Multiply B by right singular vectors
526* (CWorkspace: need N, prefer N*NRHS)
527* (RWorkspace: none)
528*
529 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
530 CALL cgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
531 $ czero, work, ldb )
532 CALL clacpy( 'G', n, nrhs, work, ldb, b, ldb )
533 ELSE IF( nrhs.GT.1 ) THEN
534 chunk = lwork / n
535 DO 20 i = 1, nrhs, chunk
536 bl = min( nrhs-i+1, chunk )
537 CALL cgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
538 $ ldb, czero, work, n )
539 CALL clacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
540 20 CONTINUE
541 ELSE IF( nrhs.EQ.1 ) THEN
542 CALL cgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
543 CALL ccopy( n, work, 1, b, 1 )
544 END IF
545*
546 ELSE IF( n.GE.mnthr .AND. lwork.GE.3*m+m*m+max( m, nrhs, n-2*m ) )
547 $ THEN
548*
549* Underdetermined case, M much less than N
550*
551* Path 2a - underdetermined, with many more columns than rows
552* and sufficient workspace for an efficient algorithm
553*
554 ldwork = m
555 IF( lwork.GE.3*m+m*lda+max( m, nrhs, n-2*m ) )
556 $ ldwork = lda
557 itau = 1
558 iwork = m + 1
559*
560* Compute A=L*Q
561* (CWorkspace: need 2*M, prefer M+M*NB)
562* (RWorkspace: none)
563*
564 CALL cgelqf( m, n, a, lda, work( itau ), work( iwork ),
565 $ lwork-iwork+1, info )
566 il = iwork
567*
568* Copy L to WORK(IL), zeroing out above it
569*
570 CALL clacpy( 'L', m, m, a, lda, work( il ), ldwork )
571 CALL claset( 'U', m-1, m-1, czero, czero, work( il+ldwork ),
572 $ ldwork )
573 ie = 1
574 itauq = il + ldwork*m
575 itaup = itauq + m
576 iwork = itaup + m
577*
578* Bidiagonalize L in WORK(IL)
579* (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
580* (RWorkspace: need M)
581*
582 CALL cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
583 $ work( itauq ), work( itaup ), work( iwork ),
584 $ lwork-iwork+1, info )
585*
586* Multiply B by transpose of left bidiagonalizing vectors of L
587* (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
588* (RWorkspace: none)
589*
590 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, work( il ), ldwork,
591 $ work( itauq ), b, ldb, work( iwork ),
592 $ lwork-iwork+1, info )
593*
594* Generate right bidiagonalizing vectors of R in WORK(IL)
595* (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
596* (RWorkspace: none)
597*
598 CALL cungbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
599 $ work( iwork ), lwork-iwork+1, info )
600 irwork = ie + m
601*
602* Perform bidiagonal QR iteration, computing right singular
603* vectors of L in WORK(IL) and multiplying B by transpose of
604* left singular vectors
605* (CWorkspace: need M*M)
606* (RWorkspace: need BDSPAC)
607*
608 CALL cbdsqr( 'U', m, m, 0, nrhs, s, rwork( ie ), work( il ),
609 $ ldwork, a, lda, b, ldb, rwork( irwork ), info )
610 IF( info.NE.0 )
611 $ GO TO 70
612*
613* Multiply B by reciprocals of singular values
614*
615 thr = max( rcond*s( 1 ), sfmin )
616 IF( rcond.LT.zero )
617 $ thr = max( eps*s( 1 ), sfmin )
618 rank = 0
619 DO 30 i = 1, m
620 IF( s( i ).GT.thr ) THEN
621 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
622 rank = rank + 1
623 ELSE
624 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
625 END IF
626 30 CONTINUE
627 iwork = il + m*ldwork
628*
629* Multiply B by right singular vectors of L in WORK(IL)
630* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
631* (RWorkspace: none)
632*
633 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
634 CALL cgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
635 $ b, ldb, czero, work( iwork ), ldb )
636 CALL clacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
637 ELSE IF( nrhs.GT.1 ) THEN
638 chunk = ( lwork-iwork+1 ) / m
639 DO 40 i = 1, nrhs, chunk
640 bl = min( nrhs-i+1, chunk )
641 CALL cgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
642 $ b( 1, i ), ldb, czero, work( iwork ), m )
643 CALL clacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
644 $ ldb )
645 40 CONTINUE
646 ELSE IF( nrhs.EQ.1 ) THEN
647 CALL cgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
648 $ 1, czero, work( iwork ), 1 )
649 CALL ccopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
650 END IF
651*
652* Zero out below first M rows of B
653*
654 CALL claset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
655 iwork = itau + m
656*
657* Multiply transpose(Q) by B
658* (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
659* (RWorkspace: none)
660*
661 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, work( itau ), b,
662 $ ldb, work( iwork ), lwork-iwork+1, info )
663*
664 ELSE
665*
666* Path 2 - remaining underdetermined cases
667*
668 ie = 1
669 itauq = 1
670 itaup = itauq + m
671 iwork = itaup + m
672*
673* Bidiagonalize A
674* (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
675* (RWorkspace: need N)
676*
677 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
678 $ work( itaup ), work( iwork ), lwork-iwork+1,
679 $ info )
680*
681* Multiply B by transpose of left bidiagonalizing vectors
682* (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
683* (RWorkspace: none)
684*
685 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda, work( itauq ),
686 $ b, ldb, work( iwork ), lwork-iwork+1, info )
687*
688* Generate right bidiagonalizing vectors in A
689* (CWorkspace: need 3*M, prefer 2*M+M*NB)
690* (RWorkspace: none)
691*
692 CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
693 $ work( iwork ), lwork-iwork+1, info )
694 irwork = ie + m
695*
696* Perform bidiagonal QR iteration,
697* computing right singular vectors of A in A and
698* multiplying B by transpose of left singular vectors
699* (CWorkspace: none)
700* (RWorkspace: need BDSPAC)
701*
702 CALL cbdsqr( 'L', m, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
703 $ 1, b, ldb, rwork( irwork ), info )
704 IF( info.NE.0 )
705 $ GO TO 70
706*
707* Multiply B by reciprocals of singular values
708*
709 thr = max( rcond*s( 1 ), sfmin )
710 IF( rcond.LT.zero )
711 $ thr = max( eps*s( 1 ), sfmin )
712 rank = 0
713 DO 50 i = 1, m
714 IF( s( i ).GT.thr ) THEN
715 CALL csrscl( nrhs, s( i ), b( i, 1 ), ldb )
716 rank = rank + 1
717 ELSE
718 CALL claset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
719 END IF
720 50 CONTINUE
721*
722* Multiply B by right singular vectors of A
723* (CWorkspace: need N, prefer N*NRHS)
724* (RWorkspace: none)
725*
726 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
727 CALL cgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
728 $ czero, work, ldb )
729 CALL clacpy( 'G', n, nrhs, work, ldb, b, ldb )
730 ELSE IF( nrhs.GT.1 ) THEN
731 chunk = lwork / n
732 DO 60 i = 1, nrhs, chunk
733 bl = min( nrhs-i+1, chunk )
734 CALL cgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
735 $ ldb, czero, work, n )
736 CALL clacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
737 60 CONTINUE
738 ELSE IF( nrhs.EQ.1 ) THEN
739 CALL cgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
740 CALL ccopy( n, work, 1, b, 1 )
741 END IF
742 END IF
743*
744* Undo scaling
745*
746 IF( iascl.EQ.1 ) THEN
747 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
748 CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
749 $ info )
750 ELSE IF( iascl.EQ.2 ) THEN
751 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
752 CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
753 $ info )
754 END IF
755 IF( ibscl.EQ.1 ) THEN
756 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
757 ELSE IF( ibscl.EQ.2 ) THEN
758 CALL clascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
759 END IF
760 70 CONTINUE
761 work( 1 ) = sroundup_lwork(maxwrk)
762 RETURN
763*
764* End of CGELSS
765*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cbdsqr
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
Definition cbdsqr.f:233
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgebrd
subroutine cgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
CGEBRD
Definition cgebrd.f:206
cgelqf
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:143
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:115
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:143
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:143
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:110
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:106
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
csrscl
subroutine csrscl(n, sa, sx, incx)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition csrscl.f:84
cungbr
subroutine cungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
CUNGBR
Definition cungbr.f:157
cunmbr
subroutine cunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMBR
Definition cunmbr.f:197
cunmlq
subroutine cunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMLQ
Definition cunmlq.f:168
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
Here is the call graph for this function:
Here is the caller graph for this function: