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

◆ zgelss()

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

ZGELSS solves overdetermined or underdetermined systems for GE matrices

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

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