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

◆ sgelss()

subroutine sgelss ( integer  m,
integer  n,
integer  nrhs,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( ldb, * )  b,
integer  ldb,
real, dimension( * )  s,
real  rcond,
integer  rank,
real, dimension( * )  work,
integer  lwork,
integer  info 
)

SGELSS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
 SGELSS computes the minimum norm solution to a real 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 REAL 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 REAL 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 elements n+1:m in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,max(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 REAL 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 >= 3*min(M,N) + max( 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]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 170 of file sgelss.f.

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