LAPACK 3.11.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, slabad, slacpy, slascl, slaset, sorgbr,
206 $ sormbr, sormlq, sormqr, srscl, xerbla
207* ..
208* .. External Functions ..
209 INTEGER ILAENV
210 REAL SLAMCH, SLANGE
211 EXTERNAL ilaenv, slamch, slange
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 ) = 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 CALL slabad( smlnum, bignum )
385*
386* Scale A if max element outside range [SMLNUM,BIGNUM]
387*
388 anrm = slange( 'M', m, n, a, lda, work )
389 iascl = 0
390 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
391*
392* Scale matrix norm up to SMLNUM
393*
394 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
395 iascl = 1
396 ELSE IF( anrm.GT.bignum ) THEN
397*
398* Scale matrix norm down to BIGNUM
399*
400 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
401 iascl = 2
402 ELSE IF( anrm.EQ.zero ) THEN
403*
404* Matrix all zero. Return zero solution.
405*
406 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
407 CALL slaset( 'F', minmn, 1, zero, zero, s, minmn )
408 rank = 0
409 GO TO 70
410 END IF
411*
412* Scale B if max element outside range [SMLNUM,BIGNUM]
413*
414 bnrm = slange( 'M', m, nrhs, b, ldb, work )
415 ibscl = 0
416 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
417*
418* Scale matrix norm up to SMLNUM
419*
420 CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
421 ibscl = 1
422 ELSE IF( bnrm.GT.bignum ) THEN
423*
424* Scale matrix norm down to BIGNUM
425*
426 CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
427 ibscl = 2
428 END IF
429*
430* Overdetermined case
431*
432 IF( m.GE.n ) THEN
433*
434* Path 1 - overdetermined or exactly determined
435*
436 mm = m
437 IF( m.GE.mnthr ) THEN
438*
439* Path 1a - overdetermined, with many more rows than columns
440*
441 mm = n
442 itau = 1
443 iwork = itau + n
444*
445* Compute A=Q*R
446* (Workspace: need 2*N, prefer N+N*NB)
447*
448 CALL sgeqrf( m, n, a, lda, work( itau ), work( iwork ),
449 $ lwork-iwork+1, info )
450*
451* Multiply B by transpose(Q)
452* (Workspace: need N+NRHS, prefer N+NRHS*NB)
453*
454 CALL sormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
455 $ ldb, work( iwork ), lwork-iwork+1, info )
456*
457* Zero out below R
458*
459 IF( n.GT.1 )
460 $ CALL slaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
461 END IF
462*
463 ie = 1
464 itauq = ie + n
465 itaup = itauq + n
466 iwork = itaup + n
467*
468* Bidiagonalize R in A
469* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
470*
471 CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
472 $ work( itaup ), work( iwork ), lwork-iwork+1,
473 $ info )
474*
475* Multiply B by transpose of left bidiagonalizing vectors of R
476* (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
477*
478 CALL sormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
479 $ b, ldb, work( iwork ), lwork-iwork+1, info )
480*
481* Generate right bidiagonalizing vectors of R in A
482* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
483*
484 CALL sorgbr( 'P', n, n, n, a, lda, work( itaup ),
485 $ work( iwork ), lwork-iwork+1, info )
486 iwork = ie + n
487*
488* Perform bidiagonal QR iteration
489* multiply B by transpose of left singular vectors
490* compute right singular vectors in A
491* (Workspace: need BDSPAC)
492*
493 CALL sbdsqr( 'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
494 $ 1, b, ldb, work( iwork ), info )
495 IF( info.NE.0 )
496 $ GO TO 70
497*
498* Multiply B by reciprocals of singular values
499*
500 thr = max( rcond*s( 1 ), sfmin )
501 IF( rcond.LT.zero )
502 $ thr = max( eps*s( 1 ), sfmin )
503 rank = 0
504 DO 10 i = 1, n
505 IF( s( i ).GT.thr ) THEN
506 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
507 rank = rank + 1
508 ELSE
509 CALL slaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
510 END IF
511 10 CONTINUE
512*
513* Multiply B by right singular vectors
514* (Workspace: need N, prefer N*NRHS)
515*
516 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
517 CALL sgemm( 'T', 'N', n, nrhs, n, one, a, lda, b, ldb, zero,
518 $ work, ldb )
519 CALL slacpy( 'G', n, nrhs, work, ldb, b, ldb )
520 ELSE IF( nrhs.GT.1 ) THEN
521 chunk = lwork / n
522 DO 20 i = 1, nrhs, chunk
523 bl = min( nrhs-i+1, chunk )
524 CALL sgemm( 'T', 'N', n, bl, n, one, a, lda, b( 1, i ),
525 $ ldb, zero, work, n )
526 CALL slacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
527 20 CONTINUE
528 ELSE
529 CALL sgemv( 'T', n, n, one, a, lda, b, 1, zero, work, 1 )
530 CALL scopy( n, work, 1, b, 1 )
531 END IF
532*
533 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
534 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
535*
536* Path 2a - underdetermined, with many more columns than rows
537* and sufficient workspace for an efficient algorithm
538*
539 ldwork = m
540 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
541 $ m*lda+m+m*nrhs ) )ldwork = lda
542 itau = 1
543 iwork = m + 1
544*
545* Compute A=L*Q
546* (Workspace: need 2*M, prefer M+M*NB)
547*
548 CALL sgelqf( m, n, a, lda, work( itau ), work( iwork ),
549 $ lwork-iwork+1, info )
550 il = iwork
551*
552* Copy L to WORK(IL), zeroing out above it
553*
554 CALL slacpy( 'L', m, m, a, lda, work( il ), ldwork )
555 CALL slaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
556 $ ldwork )
557 ie = il + ldwork*m
558 itauq = ie + m
559 itaup = itauq + m
560 iwork = itaup + m
561*
562* Bidiagonalize L in WORK(IL)
563* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
564*
565 CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
566 $ work( itauq ), work( itaup ), work( iwork ),
567 $ lwork-iwork+1, info )
568*
569* Multiply B by transpose of left bidiagonalizing vectors of L
570* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
571*
572 CALL sormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
573 $ work( itauq ), b, ldb, work( iwork ),
574 $ lwork-iwork+1, info )
575*
576* Generate right bidiagonalizing vectors of R in WORK(IL)
577* (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
578*
579 CALL sorgbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
580 $ work( iwork ), lwork-iwork+1, info )
581 iwork = ie + m
582*
583* Perform bidiagonal QR iteration,
584* computing right singular vectors of L in WORK(IL) and
585* multiplying B by transpose of left singular vectors
586* (Workspace: need M*M+M+BDSPAC)
587*
588 CALL sbdsqr( 'U', m, m, 0, nrhs, s, work( ie ), work( il ),
589 $ ldwork, a, lda, b, ldb, work( iwork ), info )
590 IF( info.NE.0 )
591 $ GO TO 70
592*
593* Multiply B by reciprocals of singular values
594*
595 thr = max( rcond*s( 1 ), sfmin )
596 IF( rcond.LT.zero )
597 $ thr = max( eps*s( 1 ), sfmin )
598 rank = 0
599 DO 30 i = 1, m
600 IF( s( i ).GT.thr ) THEN
601 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
602 rank = rank + 1
603 ELSE
604 CALL slaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
605 END IF
606 30 CONTINUE
607 iwork = ie
608*
609* Multiply B by right singular vectors of L in WORK(IL)
610* (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
611*
612 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
613 CALL sgemm( 'T', 'N', m, nrhs, m, one, work( il ), ldwork,
614 $ b, ldb, zero, work( iwork ), ldb )
615 CALL slacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
616 ELSE IF( nrhs.GT.1 ) THEN
617 chunk = ( lwork-iwork+1 ) / m
618 DO 40 i = 1, nrhs, chunk
619 bl = min( nrhs-i+1, chunk )
620 CALL sgemm( 'T', 'N', m, bl, m, one, work( il ), ldwork,
621 $ b( 1, i ), ldb, zero, work( iwork ), m )
622 CALL slacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
623 $ ldb )
624 40 CONTINUE
625 ELSE
626 CALL sgemv( 'T', m, m, one, work( il ), ldwork, b( 1, 1 ),
627 $ 1, zero, work( iwork ), 1 )
628 CALL scopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
629 END IF
630*
631* Zero out below first M rows of B
632*
633 CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
634 iwork = itau + m
635*
636* Multiply transpose(Q) by B
637* (Workspace: need M+NRHS, prefer M+NRHS*NB)
638*
639 CALL sormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
640 $ ldb, work( iwork ), lwork-iwork+1, info )
641*
642 ELSE
643*
644* Path 2 - remaining underdetermined cases
645*
646 ie = 1
647 itauq = ie + m
648 itaup = itauq + m
649 iwork = itaup + m
650*
651* Bidiagonalize A
652* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
653*
654 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
655 $ work( itaup ), work( iwork ), lwork-iwork+1,
656 $ info )
657*
658* Multiply B by transpose of left bidiagonalizing vectors
659* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
660*
661 CALL sormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
662 $ b, ldb, work( iwork ), lwork-iwork+1, info )
663*
664* Generate right bidiagonalizing vectors in A
665* (Workspace: need 4*M, prefer 3*M+M*NB)
666*
667 CALL sorgbr( 'P', m, n, m, a, lda, work( itaup ),
668 $ work( iwork ), lwork-iwork+1, info )
669 iwork = ie + m
670*
671* Perform bidiagonal QR iteration,
672* computing right singular vectors of A in A and
673* multiplying B by transpose of left singular vectors
674* (Workspace: need BDSPAC)
675*
676 CALL sbdsqr( 'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
677 $ 1, b, ldb, work( iwork ), info )
678 IF( info.NE.0 )
679 $ GO TO 70
680*
681* Multiply B by reciprocals of singular values
682*
683 thr = max( rcond*s( 1 ), sfmin )
684 IF( rcond.LT.zero )
685 $ thr = max( eps*s( 1 ), sfmin )
686 rank = 0
687 DO 50 i = 1, m
688 IF( s( i ).GT.thr ) THEN
689 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
690 rank = rank + 1
691 ELSE
692 CALL slaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
693 END IF
694 50 CONTINUE
695*
696* Multiply B by right singular vectors of A
697* (Workspace: need N, prefer N*NRHS)
698*
699 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
700 CALL sgemm( 'T', 'N', n, nrhs, m, one, a, lda, b, ldb, zero,
701 $ work, ldb )
702 CALL slacpy( 'F', n, nrhs, work, ldb, b, ldb )
703 ELSE IF( nrhs.GT.1 ) THEN
704 chunk = lwork / n
705 DO 60 i = 1, nrhs, chunk
706 bl = min( nrhs-i+1, chunk )
707 CALL sgemm( 'T', 'N', n, bl, m, one, a, lda, b( 1, i ),
708 $ ldb, zero, work, n )
709 CALL slacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
710 60 CONTINUE
711 ELSE
712 CALL sgemv( 'T', m, n, one, a, lda, b, 1, zero, work, 1 )
713 CALL scopy( n, work, 1, b, 1 )
714 END IF
715 END IF
716*
717* Undo scaling
718*
719 IF( iascl.EQ.1 ) THEN
720 CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
721 CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
722 $ info )
723 ELSE IF( iascl.EQ.2 ) THEN
724 CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
725 CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
726 $ info )
727 END IF
728 IF( ibscl.EQ.1 ) THEN
729 CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
730 ELSE IF( ibscl.EQ.2 ) THEN
731 CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
732 END IF
733*
734 70 CONTINUE
735 work( 1 ) = maxwrk
736 RETURN
737*
738* End of SGELSS
739*
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
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
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
sbdsqr
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:240
sorgbr
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
Definition: sorgbr.f:157
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
sgeqrf
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:146
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
srscl
subroutine srscl(N, SA, SX, INCX)
SRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: srscl.f:84
sormbr
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
Definition: sormbr.f:196
sormqr
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
sormlq
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:168
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
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: