LAPACK 3.11.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, slabad, slascl, slaset,
218 $ xerbla
219* ..
220* .. External Functions ..
221 INTEGER ILAENV
222 REAL CLANGE, SLAMCH
223 EXTERNAL ilaenv, clange, slamch
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, 'CGELSS', ' ', 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 CGEQRF
268 CALL cgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
269 lwork_cgeqrf = int( dum(1) )
270* Compute space needed for CUNMQR
271 CALL cunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
272 $ ldb, dum(1), -1, info )
273 lwork_cunmqr = int( dum(1) )
274 mm = n
275 maxwrk = max( maxwrk, n + n*ilaenv( 1, 'CGEQRF', ' ', m,
276 $ n, -1, -1 ) )
277 maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'CUNMQR', '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 CGEBRD
285 CALL cgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
286 $ -1, info )
287 lwork_cgebrd = int( dum(1) )
288* Compute space needed for CUNMBR
289 CALL cunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
290 $ b, ldb, dum(1), -1, info )
291 lwork_cunmbr = int( dum(1) )
292* Compute space needed for CUNGBR
293 CALL cungbr( 'P', n, n, n, a, lda, dum(1),
294 $ dum(1), -1, info )
295 lwork_cungbr = int( dum(1) )
296* Compute total workspace needed
297 maxwrk = max( maxwrk, 2*n + lwork_cgebrd )
298 maxwrk = max( maxwrk, 2*n + lwork_cunmbr )
299 maxwrk = max( maxwrk, 2*n + lwork_cungbr )
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 CGELQF
311 CALL cgelqf( m, n, a, lda, dum(1), dum(1),
312 $ -1, info )
313 lwork_cgelqf = int( dum(1) )
314* Compute space needed for CGEBRD
315 CALL cgebrd( m, m, a, lda, s, s, dum(1), dum(1),
316 $ dum(1), -1, info )
317 lwork_cgebrd = int( dum(1) )
318* Compute space needed for CUNMBR
319 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
320 $ dum(1), b, ldb, dum(1), -1, info )
321 lwork_cunmbr = int( dum(1) )
322* Compute space needed for CUNGBR
323 CALL cungbr( 'P', m, m, m, a, lda, dum(1),
324 $ dum(1), -1, info )
325 lwork_cungbr = int( dum(1) )
326* Compute space needed for CUNMLQ
327 CALL cunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
328 $ b, ldb, dum(1), -1, info )
329 lwork_cunmlq = int( dum(1) )
330* Compute total workspace needed
331 maxwrk = m + lwork_cgelqf
332 maxwrk = max( maxwrk, 3*m + m*m + lwork_cgebrd )
333 maxwrk = max( maxwrk, 3*m + m*m + lwork_cunmbr )
334 maxwrk = max( maxwrk, 3*m + m*m + lwork_cungbr )
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_cunmlq )
341 ELSE
342*
343* Path 2 - underdetermined
344*
345* Compute space needed for CGEBRD
346 CALL cgebrd( m, n, a, lda, s, s, dum(1), dum(1),
347 $ dum(1), -1, info )
348 lwork_cgebrd = int( dum(1) )
349* Compute space needed for CUNMBR
350 CALL cunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
351 $ dum(1), b, ldb, dum(1), -1, info )
352 lwork_cunmbr = int( dum(1) )
353* Compute space needed for CUNGBR
354 CALL cungbr( 'P', m, n, m, a, lda, dum(1),
355 $ dum(1), -1, info )
356 lwork_cungbr = int( dum(1) )
357 maxwrk = 2*m + lwork_cgebrd
358 maxwrk = max( maxwrk, 2*m + lwork_cunmbr )
359 maxwrk = max( maxwrk, 2*m + lwork_cungbr )
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( 'CGELSS', -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 = slamch( 'P' )
388 sfmin = slamch( 'S' )
389 smlnum = sfmin / eps
390 bignum = one / smlnum
391 CALL slabad( smlnum, bignum )
392*
393* Scale A if max element outside range [SMLNUM,BIGNUM]
394*
395 anrm = clange( '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 clascl( '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 clascl( '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 claset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
414 CALL slaset( '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 = clange( '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 clascl( '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 clascl( '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 cgeqrf( 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 cunmqr( '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 claset( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
521 rank = rank + 1
522 ELSE
523 CALL claset( '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 cgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
533 $ czero, work, ldb )
534 CALL clacpy( '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 cgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
540 $ ldb, czero, work, n )
541 CALL clacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
542 20 CONTINUE
543 ELSE
544 CALL cgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
545 CALL ccopy( 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 cgelqf( 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 clacpy( 'L', m, m, a, lda, work( il ), ldwork )
573 CALL claset( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
624 rank = rank + 1
625 ELSE
626 CALL claset( '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 cgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
637 $ b, ldb, czero, work( iwork ), ldb )
638 CALL clacpy( '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 cgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
644 $ b( 1, i ), ldb, czero, work( iwork ), m )
645 CALL clacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
646 $ ldb )
647 40 CONTINUE
648 ELSE
649 CALL cgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
650 $ 1, czero, work( iwork ), 1 )
651 CALL ccopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
652 END IF
653*
654* Zero out below first M rows of B
655*
656 CALL claset( '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 cunmlq( '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 cgebrd( 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 cunmbr( '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 cungbr( '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 cbdsqr( '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 csrscl( nrhs, s( i ), b( i, 1 ), ldb )
718 rank = rank + 1
719 ELSE
720 CALL claset( '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 cgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
730 $ czero, work, ldb )
731 CALL clacpy( '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 cgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
737 $ ldb, czero, work, n )
738 CALL clacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
739 60 CONTINUE
740 ELSE
741 CALL cgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
742 CALL ccopy( n, work, 1, b, 1 )
743 END IF
744 END IF
745*
746* Undo scaling
747*
748 IF( iascl.EQ.1 ) THEN
749 CALL clascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
750 CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
751 $ info )
752 ELSE IF( iascl.EQ.2 ) THEN
753 CALL clascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
754 CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
755 $ info )
756 END IF
757 IF( ibscl.EQ.1 ) THEN
758 CALL clascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
759 ELSE IF( ibscl.EQ.2 ) THEN
760 CALL clascl( '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 CGELSS
767*
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
cungbr
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
Definition: cungbr.f:157
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
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
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
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
csrscl
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: csrscl.f:84
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
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
cunmbr
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:197
cbdsqr
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
Definition: cbdsqr.f:222
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
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: