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

◆ dgelss()

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

DGELSS solves overdetermined or underdetermined systems for GE matrices

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

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