LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zgelss.f
Go to the documentation of this file.
1*> \brief <b> ZGELSS solves overdetermined or underdetermined systems for GE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelss.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelss.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelss.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
22* WORK, LWORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26* DOUBLE PRECISION RCOND
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION RWORK( * ), S( * )
30* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> ZGELSS computes the minimum norm solution to a complex linear
40*> least squares problem:
41*>
42*> Minimize 2-norm(| b - A*x |).
43*>
44*> using the singular value decomposition (SVD) of A. A is an M-by-N
45*> matrix which may be rank-deficient.
46*>
47*> Several right hand side vectors b and solution vectors x can be
48*> handled in a single call; they are stored as the columns of the
49*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
50*> X.
51*>
52*> The effective rank of A is determined by treating as zero those
53*> singular values which are less than RCOND times the largest singular
54*> value.
55*> \endverbatim
56*
57* Arguments:
58* ==========
59*
60*> \param[in] M
61*> \verbatim
62*> M is INTEGER
63*> The number of rows of the matrix A. M >= 0.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*> N is INTEGER
69*> The number of columns of the matrix A. N >= 0.
70*> \endverbatim
71*>
72*> \param[in] NRHS
73*> \verbatim
74*> NRHS is INTEGER
75*> The number of right hand sides, i.e., the number of columns
76*> of the matrices B and X. NRHS >= 0.
77*> \endverbatim
78*>
79*> \param[in,out] A
80*> \verbatim
81*> A is COMPLEX*16 array, dimension (LDA,N)
82*> On entry, the M-by-N matrix A.
83*> On exit, the first min(m,n) rows of A are overwritten with
84*> its right singular vectors, stored rowwise.
85*> \endverbatim
86*>
87*> \param[in] LDA
88*> \verbatim
89*> LDA is INTEGER
90*> The leading dimension of the array A. LDA >= max(1,M).
91*> \endverbatim
92*>
93*> \param[in,out] B
94*> \verbatim
95*> B is COMPLEX*16 array, dimension (LDB,NRHS)
96*> On entry, the M-by-NRHS right hand side matrix B.
97*> On exit, B is overwritten by the N-by-NRHS solution matrix X.
98*> If m >= n and RANK = n, the residual sum-of-squares for
99*> the solution in the i-th column is given by the sum of
100*> squares of the modulus of elements n+1:m in that column.
101*> \endverbatim
102*>
103*> \param[in] LDB
104*> \verbatim
105*> LDB is INTEGER
106*> The leading dimension of the array B. LDB >= max(1,M,N).
107*> \endverbatim
108*>
109*> \param[out] S
110*> \verbatim
111*> S is DOUBLE PRECISION array, dimension (min(M,N))
112*> The singular values of A in decreasing order.
113*> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
114*> \endverbatim
115*>
116*> \param[in] RCOND
117*> \verbatim
118*> RCOND is DOUBLE PRECISION
119*> RCOND is used to determine the effective rank of A.
120*> Singular values S(i) <= RCOND*S(1) are treated as zero.
121*> If RCOND < 0, machine precision is used instead.
122*> \endverbatim
123*>
124*> \param[out] RANK
125*> \verbatim
126*> RANK is INTEGER
127*> The effective rank of A, i.e., the number of singular values
128*> which are greater than RCOND*S(1).
129*> \endverbatim
130*>
131*> \param[out] WORK
132*> \verbatim
133*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
134*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
135*> \endverbatim
136*>
137*> \param[in] LWORK
138*> \verbatim
139*> LWORK is INTEGER
140*> The dimension of the array WORK. LWORK >= 1, and also:
141*> LWORK >= 2*min(M,N) + max(M,N,NRHS)
142*> For good performance, LWORK should generally be larger.
143*>
144*> If LWORK = -1, then a workspace query is assumed; the routine
145*> only calculates the optimal size of the WORK array, returns
146*> this value as the first entry of the WORK array, and no error
147*> message related to LWORK is issued by XERBLA.
148*> \endverbatim
149*>
150*> \param[out] RWORK
151*> \verbatim
152*> RWORK is DOUBLE PRECISION array, dimension (5*min(M,N))
153*> \endverbatim
154*>
155*> \param[out] INFO
156*> \verbatim
157*> INFO is INTEGER
158*> = 0: successful exit
159*> < 0: if INFO = -i, the i-th argument had an illegal value.
160*> > 0: the algorithm for computing the SVD failed to converge;
161*> if INFO = i, i off-diagonal elements of an intermediate
162*> bidiagonal form did not converge to zero.
163*> \endverbatim
164*
165* Authors:
166* ========
167*
168*> \author Univ. of Tennessee
169*> \author Univ. of California Berkeley
170*> \author Univ. of Colorado Denver
171*> \author NAG Ltd.
172*
173*> \ingroup complex16GEsolve
174*
175* =====================================================================
176 SUBROUTINE zgelss( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
177 \$ WORK, LWORK, RWORK, INFO )
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 dlabad, dlascl, dlaset, xerbla, zbdsqr, zcopy,
218 \$ zunmqr
219* ..
220* .. External Functions ..
221 INTEGER ILAENV
222 DOUBLE PRECISION DLAMCH, ZLANGE
223 EXTERNAL ilaenv, dlamch, zlange
224* ..
225* .. Intrinsic Functions ..
226 INTRINSIC max, min
227* ..
228* .. Executable Statements ..
229*
230* Test the input arguments
231*
232 info = 0
233 minmn = min( m, n )
234 maxmn = max( m, n )
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( nrhs.LT.0 ) THEN
241 info = -3
242 ELSE IF( lda.LT.max( 1, m ) ) THEN
243 info = -5
244 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
245 info = -7
246 END IF
247*
248* Compute workspace
249* (Note: Comments in the code beginning "Workspace:" describe the
250* minimal amount of workspace needed at that point in the code,
251* as well as the preferred amount for good performance.
252* CWorkspace refers to complex workspace, and RWorkspace refers
253* to real workspace. NB refers to the optimal block size for the
254* immediately following subroutine, as returned by ILAENV.)
255*
256 IF( info.EQ.0 ) THEN
257 minwrk = 1
258 maxwrk = 1
259 IF( minmn.GT.0 ) THEN
260 mm = m
261 mnthr = ilaenv( 6, 'ZGELSS', ' ', m, n, nrhs, -1 )
262 IF( m.GE.n .AND. m.GE.mnthr ) THEN
263*
264* Path 1a - overdetermined, with many more rows than
265* columns
266*
267* Compute space needed for ZGEQRF
268 CALL zgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
269 lwork_zgeqrf = int( dum(1) )
270* Compute space needed for ZUNMQR
271 CALL zunmqr( 'L', 'C', m, nrhs, n, a, lda, dum(1), b,
272 \$ ldb, dum(1), -1, info )
273 lwork_zunmqr = int( dum(1) )
274 mm = n
275 maxwrk = max( maxwrk, n + n*ilaenv( 1, 'ZGEQRF', ' ', m,
276 \$ n, -1, -1 ) )
277 maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'ZUNMQR', 'LC',
278 \$ m, nrhs, n, -1 ) )
279 END IF
280 IF( m.GE.n ) THEN
281*
282* Path 1 - overdetermined or exactly determined
283*
284* Compute space needed for ZGEBRD
285 CALL zgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
286 \$ -1, info )
287 lwork_zgebrd = int( dum(1) )
288* Compute space needed for ZUNMBR
289 CALL zunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, dum(1),
290 \$ b, ldb, dum(1), -1, info )
291 lwork_zunmbr = int( dum(1) )
292* Compute space needed for ZUNGBR
293 CALL zungbr( 'P', n, n, n, a, lda, dum(1),
294 \$ dum(1), -1, info )
295 lwork_zungbr = int( dum(1) )
296* Compute total workspace needed
297 maxwrk = max( maxwrk, 2*n + lwork_zgebrd )
298 maxwrk = max( maxwrk, 2*n + lwork_zunmbr )
299 maxwrk = max( maxwrk, 2*n + lwork_zungbr )
300 maxwrk = max( maxwrk, n*nrhs )
301 minwrk = 2*n + max( nrhs, m )
302 END IF
303 IF( n.GT.m ) THEN
304 minwrk = 2*m + max( nrhs, n )
305 IF( n.GE.mnthr ) THEN
306*
307* Path 2a - underdetermined, with many more columns
308* than rows
309*
310* Compute space needed for ZGELQF
311 CALL zgelqf( m, n, a, lda, dum(1), dum(1),
312 \$ -1, info )
313 lwork_zgelqf = int( dum(1) )
314* Compute space needed for ZGEBRD
315 CALL zgebrd( m, m, a, lda, s, s, dum(1), dum(1),
316 \$ dum(1), -1, info )
317 lwork_zgebrd = int( dum(1) )
318* Compute space needed for ZUNMBR
319 CALL zunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda,
320 \$ dum(1), b, ldb, dum(1), -1, info )
321 lwork_zunmbr = int( dum(1) )
322* Compute space needed for ZUNGBR
323 CALL zungbr( 'P', m, m, m, a, lda, dum(1),
324 \$ dum(1), -1, info )
325 lwork_zungbr = int( dum(1) )
326* Compute space needed for ZUNMLQ
327 CALL zunmlq( 'L', 'C', n, nrhs, m, a, lda, dum(1),
328 \$ b, ldb, dum(1), -1, info )
329 lwork_zunmlq = int( dum(1) )
330* Compute total workspace needed
331 maxwrk = m + lwork_zgelqf
332 maxwrk = max( maxwrk, 3*m + m*m + lwork_zgebrd )
333 maxwrk = max( maxwrk, 3*m + m*m + lwork_zunmbr )
334 maxwrk = max( maxwrk, 3*m + m*m + lwork_zungbr )
335 IF( nrhs.GT.1 ) THEN
336 maxwrk = max( maxwrk, m*m + m + m*nrhs )
337 ELSE
338 maxwrk = max( maxwrk, m*m + 2*m )
339 END IF
340 maxwrk = max( maxwrk, m + lwork_zunmlq )
341 ELSE
342*
343* Path 2 - underdetermined
344*
345* Compute space needed for ZGEBRD
346 CALL zgebrd( m, n, a, lda, s, s, dum(1), dum(1),
347 \$ dum(1), -1, info )
348 lwork_zgebrd = int( dum(1) )
349* Compute space needed for ZUNMBR
350 CALL zunmbr( 'Q', 'L', 'C', m, nrhs, m, a, lda,
351 \$ dum(1), b, ldb, dum(1), -1, info )
352 lwork_zunmbr = int( dum(1) )
353* Compute space needed for ZUNGBR
354 CALL zungbr( 'P', m, n, m, a, lda, dum(1),
355 \$ dum(1), -1, info )
356 lwork_zungbr = int( dum(1) )
357 maxwrk = 2*m + lwork_zgebrd
358 maxwrk = max( maxwrk, 2*m + lwork_zunmbr )
359 maxwrk = max( maxwrk, 2*m + lwork_zungbr )
360 maxwrk = max( maxwrk, n*nrhs )
361 END IF
362 END IF
363 maxwrk = max( minwrk, maxwrk )
364 END IF
365 work( 1 ) = maxwrk
366*
367 IF( lwork.LT.minwrk .AND. .NOT.lquery )
368 \$ info = -12
369 END IF
370*
371 IF( info.NE.0 ) THEN
372 CALL xerbla( 'ZGELSS', -info )
373 RETURN
374 ELSE IF( lquery ) THEN
375 RETURN
376 END IF
377*
378* Quick return if possible
379*
380 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
381 rank = 0
382 RETURN
383 END IF
384*
385* Get machine parameters
386*
387 eps = dlamch( 'P' )
388 sfmin = dlamch( 'S' )
389 smlnum = sfmin / eps
390 bignum = one / smlnum
391 CALL dlabad( smlnum, bignum )
392*
393* Scale A if max element outside range [SMLNUM,BIGNUM]
394*
395 anrm = zlange( 'M', m, n, a, lda, rwork )
396 iascl = 0
397 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
398*
399* Scale matrix norm up to SMLNUM
400*
401 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
402 iascl = 1
403 ELSE IF( anrm.GT.bignum ) THEN
404*
405* Scale matrix norm down to BIGNUM
406*
407 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
408 iascl = 2
409 ELSE IF( anrm.EQ.zero ) THEN
410*
411* Matrix all zero. Return zero solution.
412*
413 CALL zlaset( 'F', max( m, n ), nrhs, czero, czero, b, ldb )
414 CALL dlaset( 'F', minmn, 1, zero, zero, s, minmn )
415 rank = 0
416 GO TO 70
417 END IF
418*
419* Scale B if max element outside range [SMLNUM,BIGNUM]
420*
421 bnrm = zlange( 'M', m, nrhs, b, ldb, rwork )
422 ibscl = 0
423 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
424*
425* Scale matrix norm up to SMLNUM
426*
427 CALL zlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
428 ibscl = 1
429 ELSE IF( bnrm.GT.bignum ) THEN
430*
431* Scale matrix norm down to BIGNUM
432*
433 CALL zlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
434 ibscl = 2
435 END IF
436*
437* Overdetermined case
438*
439 IF( m.GE.n ) THEN
440*
441* Path 1 - overdetermined or exactly determined
442*
443 mm = m
444 IF( m.GE.mnthr ) THEN
445*
446* Path 1a - overdetermined, with many more rows than columns
447*
448 mm = n
449 itau = 1
450 iwork = itau + n
451*
452* Compute A=Q*R
453* (CWorkspace: need 2*N, prefer N+N*NB)
454* (RWorkspace: none)
455*
456 CALL zgeqrf( m, n, a, lda, work( itau ), work( iwork ),
457 \$ lwork-iwork+1, info )
458*
459* Multiply B by transpose(Q)
460* (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
461* (RWorkspace: none)
462*
463 CALL zunmqr( 'L', 'C', m, nrhs, n, a, lda, work( itau ), b,
464 \$ ldb, work( iwork ), lwork-iwork+1, info )
465*
466* Zero out below R
467*
468 IF( n.GT.1 )
469 \$ CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
470 \$ lda )
471 END IF
472*
473 ie = 1
474 itauq = 1
475 itaup = itauq + n
476 iwork = itaup + n
477*
478* Bidiagonalize R in A
479* (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
480* (RWorkspace: need N)
481*
482 CALL zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
483 \$ work( itaup ), work( iwork ), lwork-iwork+1,
484 \$ info )
485*
486* Multiply B by transpose of left bidiagonalizing vectors of R
487* (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
488* (RWorkspace: none)
489*
490 CALL zunmbr( 'Q', 'L', 'C', mm, nrhs, n, a, lda, work( itauq ),
491 \$ b, ldb, work( iwork ), lwork-iwork+1, info )
492*
493* Generate right bidiagonalizing vectors of R in A
494* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
495* (RWorkspace: none)
496*
497 CALL zungbr( 'P', n, n, n, a, lda, work( itaup ),
498 \$ work( iwork ), lwork-iwork+1, info )
499 irwork = ie + n
500*
501* Perform bidiagonal QR iteration
502* multiply B by transpose of left singular vectors
503* compute right singular vectors in A
504* (CWorkspace: none)
505* (RWorkspace: need BDSPAC)
506*
507 CALL zbdsqr( 'U', n, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
508 \$ 1, b, ldb, rwork( irwork ), info )
509 IF( info.NE.0 )
510 \$ GO TO 70
511*
512* Multiply B by reciprocals of singular values
513*
514 thr = max( rcond*s( 1 ), sfmin )
515 IF( rcond.LT.zero )
516 \$ thr = max( eps*s( 1 ), sfmin )
517 rank = 0
518 DO 10 i = 1, n
519 IF( s( i ).GT.thr ) THEN
520 CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
521 rank = rank + 1
522 ELSE
523 CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
524 END IF
525 10 CONTINUE
526*
527* Multiply B by right singular vectors
528* (CWorkspace: need N, prefer N*NRHS)
529* (RWorkspace: none)
530*
531 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
532 CALL zgemm( 'C', 'N', n, nrhs, n, cone, a, lda, b, ldb,
533 \$ czero, work, ldb )
534 CALL zlacpy( 'G', n, nrhs, work, ldb, b, ldb )
535 ELSE IF( nrhs.GT.1 ) THEN
536 chunk = lwork / n
537 DO 20 i = 1, nrhs, chunk
538 bl = min( nrhs-i+1, chunk )
539 CALL zgemm( 'C', 'N', n, bl, n, cone, a, lda, b( 1, i ),
540 \$ ldb, czero, work, n )
541 CALL zlacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
542 20 CONTINUE
543 ELSE
544 CALL zgemv( 'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
545 CALL zcopy( n, work, 1, b, 1 )
546 END IF
547*
548 ELSE IF( n.GE.mnthr .AND. lwork.GE.3*m+m*m+max( m, nrhs, n-2*m ) )
549 \$ THEN
550*
551* Underdetermined case, M much less than N
552*
553* Path 2a - underdetermined, with many more columns than rows
554* and sufficient workspace for an efficient algorithm
555*
556 ldwork = m
557 IF( lwork.GE.3*m+m*lda+max( m, nrhs, n-2*m ) )
558 \$ ldwork = lda
559 itau = 1
560 iwork = m + 1
561*
562* Compute A=L*Q
563* (CWorkspace: need 2*M, prefer M+M*NB)
564* (RWorkspace: none)
565*
566 CALL zgelqf( m, n, a, lda, work( itau ), work( iwork ),
567 \$ lwork-iwork+1, info )
568 il = iwork
569*
570* Copy L to WORK(IL), zeroing out above it
571*
572 CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwork )
573 CALL zlaset( 'U', m-1, m-1, czero, czero, work( il+ldwork ),
574 \$ ldwork )
575 ie = 1
576 itauq = il + ldwork*m
577 itaup = itauq + m
578 iwork = itaup + m
579*
580* Bidiagonalize L in WORK(IL)
581* (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
582* (RWorkspace: need M)
583*
584 CALL zgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
585 \$ work( itauq ), work( itaup ), work( iwork ),
586 \$ lwork-iwork+1, info )
587*
588* Multiply B by transpose of left bidiagonalizing vectors of L
589* (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
590* (RWorkspace: none)
591*
592 CALL zunmbr( 'Q', 'L', 'C', m, nrhs, m, work( il ), ldwork,
593 \$ work( itauq ), b, ldb, work( iwork ),
594 \$ lwork-iwork+1, info )
595*
596* Generate right bidiagonalizing vectors of R in WORK(IL)
597* (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
598* (RWorkspace: none)
599*
600 CALL zungbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
601 \$ work( iwork ), lwork-iwork+1, info )
602 irwork = ie + m
603*
604* Perform bidiagonal QR iteration, computing right singular
605* vectors of L in WORK(IL) and multiplying B by transpose of
606* left singular vectors
607* (CWorkspace: need M*M)
608* (RWorkspace: need BDSPAC)
609*
610 CALL zbdsqr( 'U', m, m, 0, nrhs, s, rwork( ie ), work( il ),
611 \$ ldwork, a, lda, b, ldb, rwork( irwork ), info )
612 IF( info.NE.0 )
613 \$ GO TO 70
614*
615* Multiply B by reciprocals of singular values
616*
617 thr = max( rcond*s( 1 ), sfmin )
618 IF( rcond.LT.zero )
619 \$ thr = max( eps*s( 1 ), sfmin )
620 rank = 0
621 DO 30 i = 1, m
622 IF( s( i ).GT.thr ) THEN
623 CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
624 rank = rank + 1
625 ELSE
626 CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
627 END IF
628 30 CONTINUE
629 iwork = il + m*ldwork
630*
631* Multiply B by right singular vectors of L in WORK(IL)
632* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
633* (RWorkspace: none)
634*
635 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
636 CALL zgemm( 'C', 'N', m, nrhs, m, cone, work( il ), ldwork,
637 \$ b, ldb, czero, work( iwork ), ldb )
638 CALL zlacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
639 ELSE IF( nrhs.GT.1 ) THEN
640 chunk = ( lwork-iwork+1 ) / m
641 DO 40 i = 1, nrhs, chunk
642 bl = min( nrhs-i+1, chunk )
643 CALL zgemm( 'C', 'N', m, bl, m, cone, work( il ), ldwork,
644 \$ b( 1, i ), ldb, czero, work( iwork ), m )
645 CALL zlacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
646 \$ ldb )
647 40 CONTINUE
648 ELSE
649 CALL zgemv( 'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
650 \$ 1, czero, work( iwork ), 1 )
651 CALL zcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
652 END IF
653*
654* Zero out below first M rows of B
655*
656 CALL zlaset( 'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
657 iwork = itau + m
658*
659* Multiply transpose(Q) by B
660* (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
661* (RWorkspace: none)
662*
663 CALL zunmlq( 'L', 'C', n, nrhs, m, a, lda, work( itau ), b,
664 \$ ldb, work( iwork ), lwork-iwork+1, info )
665*
666 ELSE
667*
668* Path 2 - remaining underdetermined cases
669*
670 ie = 1
671 itauq = 1
672 itaup = itauq + m
673 iwork = itaup + m
674*
675* Bidiagonalize A
676* (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
677* (RWorkspace: need N)
678*
679 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
680 \$ work( itaup ), work( iwork ), lwork-iwork+1,
681 \$ info )
682*
683* Multiply B by transpose of left bidiagonalizing vectors
684* (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
685* (RWorkspace: none)
686*
687 CALL zunmbr( 'Q', 'L', 'C', m, nrhs, n, a, lda, work( itauq ),
688 \$ b, ldb, work( iwork ), lwork-iwork+1, info )
689*
690* Generate right bidiagonalizing vectors in A
691* (CWorkspace: need 3*M, prefer 2*M+M*NB)
692* (RWorkspace: none)
693*
694 CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
695 \$ work( iwork ), lwork-iwork+1, info )
696 irwork = ie + m
697*
698* Perform bidiagonal QR iteration,
699* computing right singular vectors of A in A and
700* multiplying B by transpose of left singular vectors
701* (CWorkspace: none)
702* (RWorkspace: need BDSPAC)
703*
704 CALL zbdsqr( 'L', m, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
705 \$ 1, b, ldb, rwork( irwork ), info )
706 IF( info.NE.0 )
707 \$ GO TO 70
708*
709* Multiply B by reciprocals of singular values
710*
711 thr = max( rcond*s( 1 ), sfmin )
712 IF( rcond.LT.zero )
713 \$ thr = max( eps*s( 1 ), sfmin )
714 rank = 0
715 DO 50 i = 1, m
716 IF( s( i ).GT.thr ) THEN
717 CALL zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
718 rank = rank + 1
719 ELSE
720 CALL zlaset( 'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
721 END IF
722 50 CONTINUE
723*
724* Multiply B by right singular vectors of A
725* (CWorkspace: need N, prefer N*NRHS)
726* (RWorkspace: none)
727*
728 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
729 CALL zgemm( 'C', 'N', n, nrhs, m, cone, a, lda, b, ldb,
730 \$ czero, work, ldb )
731 CALL zlacpy( 'G', n, nrhs, work, ldb, b, ldb )
732 ELSE IF( nrhs.GT.1 ) THEN
733 chunk = lwork / n
734 DO 60 i = 1, nrhs, chunk
735 bl = min( nrhs-i+1, chunk )
736 CALL zgemm( 'C', 'N', n, bl, m, cone, a, lda, b( 1, i ),
737 \$ ldb, czero, work, n )
738 CALL zlacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
739 60 CONTINUE
740 ELSE
741 CALL zgemv( 'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
742 CALL zcopy( n, work, 1, b, 1 )
743 END IF
744 END IF
745*
746* Undo scaling
747*
748 IF( iascl.EQ.1 ) THEN
749 CALL zlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
750 CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
751 \$ info )
752 ELSE IF( iascl.EQ.2 ) THEN
753 CALL zlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
754 CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
755 \$ info )
756 END IF
757 IF( ibscl.EQ.1 ) THEN
758 CALL zlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
759 ELSE IF( ibscl.EQ.2 ) THEN
760 CALL zlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
761 END IF
762 70 CONTINUE
763 work( 1 ) = maxwrk
764 RETURN
765*
766* End of ZGELSS
767*
768 END
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
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:157
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:143
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:205
subroutine zgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSS solves overdetermined or underdetermined systems for GE matrices
Definition: zgelss.f:178
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
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
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
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:84
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:167
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
Definition: zbdsqr.f:222
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:196
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:152