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
9: *> Download ZGELSS + dependencies
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: *> \date June 2016
174: *
175: *> \ingroup complex16GEsolve
176: *
177: * =====================================================================
178: SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
179: $ WORK, LWORK, RWORK, INFO )
180: *
181: * -- LAPACK driver routine (version 3.6.1) --
182: * -- LAPACK is a software package provided by Univ. of Tennessee, --
183: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184: * June 2016
185: *
186: * .. Scalar Arguments ..
187: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
188: DOUBLE PRECISION RCOND
189: * ..
190: * .. Array Arguments ..
191: DOUBLE PRECISION RWORK( * ), S( * )
192: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
193: * ..
194: *
195: * =====================================================================
196: *
197: * .. Parameters ..
198: DOUBLE PRECISION ZERO, ONE
199: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
200: COMPLEX*16 CZERO, CONE
201: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
202: $ CONE = ( 1.0D+0, 0.0D+0 ) )
203: * ..
204: * .. Local Scalars ..
205: LOGICAL LQUERY
206: INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
207: $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
208: $ MAXWRK, MINMN, MINWRK, MM, MNTHR
209: INTEGER LWORK_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD,
210: $ LWORK_ZUNMBR, LWORK_ZUNGBR, LWORK_ZUNMLQ,
211: $ LWORK_ZGELQF
212: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
213: * ..
214: * .. Local Arrays ..
215: COMPLEX*16 DUM( 1 )
216: * ..
217: * .. External Subroutines ..
218: EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY,
219: $ ZDRSCL, ZGEBRD, ZGELQF, ZGEMM, ZGEMV, ZGEQRF,
220: $ ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNMBR, ZUNMLQ,
221: $ ZUNMQR
222: * ..
223: * .. External Functions ..
224: INTEGER ILAENV
225: DOUBLE PRECISION DLAMCH, ZLANGE
226: EXTERNAL ILAENV, DLAMCH, ZLANGE
227: * ..
228: * .. Intrinsic Functions ..
229: INTRINSIC MAX, MIN
230: * ..
231: * .. Executable Statements ..
232: *
233: * Test the input arguments
234: *
235: INFO = 0
236: MINMN = MIN( M, N )
237: MAXMN = MAX( M, N )
238: LQUERY = ( LWORK.EQ.-1 )
239: IF( M.LT.0 ) THEN
240: INFO = -1
241: ELSE IF( N.LT.0 ) THEN
242: INFO = -2
243: ELSE IF( NRHS.LT.0 ) THEN
244: INFO = -3
245: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
246: INFO = -5
247: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
248: INFO = -7
249: END IF
250: *
251: * Compute workspace
252: * (Note: Comments in the code beginning "Workspace:" describe the
253: * minimal amount of workspace needed at that point in the code,
254: * as well as the preferred amount for good performance.
255: * CWorkspace refers to complex workspace, and RWorkspace refers
256: * to real workspace. NB refers to the optimal block size for the
257: * immediately following subroutine, as returned by ILAENV.)
258: *
259: IF( INFO.EQ.0 ) THEN
260: MINWRK = 1
261: MAXWRK = 1
262: IF( MINMN.GT.0 ) THEN
263: MM = M
264: MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
265: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
266: *
267: * Path 1a - overdetermined, with many more rows than
268: * columns
269: *
270: * Compute space needed for ZGEQRF
271: CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
272: LWORK_ZGEQRF=DUM(1)
273: * Compute space needed for ZUNMQR
274: CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B,
275: $ LDB, DUM(1), -1, INFO )
276: LWORK_ZUNMQR=DUM(1)
277: MM = N
278: MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
279: $ N, -1, -1 ) )
280: MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
281: $ M, NRHS, N, -1 ) )
282: END IF
283: IF( M.GE.N ) THEN
284: *
285: * Path 1 - overdetermined or exactly determined
286: *
287: * Compute space needed for ZGEBRD
288: CALL ZGEBRD( MM, N, A, LDA, S, S, DUM(1), DUM(1), DUM(1),
289: $ -1, INFO )
290: LWORK_ZGEBRD=DUM(1)
291: * Compute space needed for ZUNMBR
292: CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1),
293: $ B, LDB, DUM(1), -1, INFO )
294: LWORK_ZUNMBR=DUM(1)
295: * Compute space needed for ZUNGBR
296: CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1),
297: $ DUM(1), -1, INFO )
298: LWORK_ZUNGBR=DUM(1)
299: * Compute total workspace needed
300: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD )
301: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR )
302: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR )
303: MAXWRK = MAX( MAXWRK, N*NRHS )
304: MINWRK = 2*N + MAX( NRHS, M )
305: END IF
306: IF( N.GT.M ) THEN
307: MINWRK = 2*M + MAX( NRHS, N )
308: IF( N.GE.MNTHR ) THEN
309: *
310: * Path 2a - underdetermined, with many more columns
311: * than rows
312: *
313: * Compute space needed for ZGELQF
314: CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1),
315: $ -1, INFO )
316: LWORK_ZGELQF=DUM(1)
317: * Compute space needed for ZGEBRD
318: CALL ZGEBRD( M, M, A, LDA, S, S, DUM(1), DUM(1),
319: $ DUM(1), -1, INFO )
320: LWORK_ZGEBRD=DUM(1)
321: * Compute space needed for ZUNMBR
322: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA,
323: $ DUM(1), B, LDB, DUM(1), -1, INFO )
324: LWORK_ZUNMBR=DUM(1)
325: * Compute space needed for ZUNGBR
326: CALL ZUNGBR( 'P', M, M, M, A, LDA, DUM(1),
327: $ DUM(1), -1, INFO )
328: LWORK_ZUNGBR=DUM(1)
329: * Compute space needed for ZUNMLQ
330: CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1),
331: $ B, LDB, DUM(1), -1, INFO )
332: LWORK_ZUNMLQ=DUM(1)
333: * Compute total workspace needed
334: MAXWRK = M + LWORK_ZGELQF
335: MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZGEBRD )
336: MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNMBR )
337: MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNGBR )
338: IF( NRHS.GT.1 ) THEN
339: MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
340: ELSE
341: MAXWRK = MAX( MAXWRK, M*M + 2*M )
342: END IF
343: MAXWRK = MAX( MAXWRK, M + LWORK_ZUNMLQ )
344: ELSE
345: *
346: * Path 2 - underdetermined
347: *
348: * Compute space needed for ZGEBRD
349: CALL ZGEBRD( M, N, A, LDA, S, S, DUM(1), DUM(1),
350: $ DUM(1), -1, INFO )
351: LWORK_ZGEBRD=DUM(1)
352: * Compute space needed for ZUNMBR
353: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA,
354: $ DUM(1), B, LDB, DUM(1), -1, INFO )
355: LWORK_ZUNMBR=DUM(1)
356: * Compute space needed for ZUNGBR
357: CALL ZUNGBR( 'P', M, N, M, A, LDA, DUM(1),
358: $ DUM(1), -1, INFO )
359: LWORK_ZUNGBR=DUM(1)
360: MAXWRK = 2*M + LWORK_ZGEBRD
361: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR )
362: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR )
363: MAXWRK = MAX( MAXWRK, N*NRHS )
364: END IF
365: END IF
366: MAXWRK = MAX( MINWRK, MAXWRK )
367: END IF
368: WORK( 1 ) = MAXWRK
369: *
370: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
371: $ INFO = -12
372: END IF
373: *
374: IF( INFO.NE.0 ) THEN
375: CALL XERBLA( 'ZGELSS', -INFO )
376: RETURN
377: ELSE IF( LQUERY ) THEN
378: RETURN
379: END IF
380: *
381: * Quick return if possible
382: *
383: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
384: RANK = 0
385: RETURN
386: END IF
387: *
388: * Get machine parameters
389: *
390: EPS = DLAMCH( 'P' )
391: SFMIN = DLAMCH( 'S' )
392: SMLNUM = SFMIN / EPS
393: BIGNUM = ONE / SMLNUM
394: CALL DLABAD( SMLNUM, BIGNUM )
395: *
396: * Scale A if max element outside range [SMLNUM,BIGNUM]
397: *
398: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
399: IASCL = 0
400: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
401: *
402: * Scale matrix norm up to SMLNUM
403: *
404: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
405: IASCL = 1
406: ELSE IF( ANRM.GT.BIGNUM ) THEN
407: *
408: * Scale matrix norm down to BIGNUM
409: *
410: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
411: IASCL = 2
412: ELSE IF( ANRM.EQ.ZERO ) THEN
413: *
414: * Matrix all zero. Return zero solution.
415: *
416: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
417: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
418: RANK = 0
419: GO TO 70
420: END IF
421: *
422: * Scale B if max element outside range [SMLNUM,BIGNUM]
423: *
424: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
425: IBSCL = 0
426: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
427: *
428: * Scale matrix norm up to SMLNUM
429: *
430: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
431: IBSCL = 1
432: ELSE IF( BNRM.GT.BIGNUM ) THEN
433: *
434: * Scale matrix norm down to BIGNUM
435: *
436: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
437: IBSCL = 2
438: END IF
439: *
440: * Overdetermined case
441: *
442: IF( M.GE.N ) THEN
443: *
444: * Path 1 - overdetermined or exactly determined
445: *
446: MM = M
447: IF( M.GE.MNTHR ) THEN
448: *
449: * Path 1a - overdetermined, with many more rows than columns
450: *
451: MM = N
452: ITAU = 1
453: IWORK = ITAU + N
454: *
455: * Compute A=Q*R
456: * (CWorkspace: need 2*N, prefer N+N*NB)
457: * (RWorkspace: none)
458: *
459: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
460: $ LWORK-IWORK+1, INFO )
461: *
462: * Multiply B by transpose(Q)
463: * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
464: * (RWorkspace: none)
465: *
466: CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
467: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
468: *
469: * Zero out below R
470: *
471: IF( N.GT.1 )
472: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
473: $ LDA )
474: END IF
475: *
476: IE = 1
477: ITAUQ = 1
478: ITAUP = ITAUQ + N
479: IWORK = ITAUP + N
480: *
481: * Bidiagonalize R in A
482: * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
483: * (RWorkspace: need N)
484: *
485: CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
486: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
487: $ INFO )
488: *
489: * Multiply B by transpose of left bidiagonalizing vectors of R
490: * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
491: * (RWorkspace: none)
492: *
493: CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
494: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
495: *
496: * Generate right bidiagonalizing vectors of R in A
497: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
498: * (RWorkspace: none)
499: *
500: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
501: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
502: IRWORK = IE + N
503: *
504: * Perform bidiagonal QR iteration
505: * multiply B by transpose of left singular vectors
506: * compute right singular vectors in A
507: * (CWorkspace: none)
508: * (RWorkspace: need BDSPAC)
509: *
510: CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
511: $ 1, B, LDB, RWORK( IRWORK ), INFO )
512: IF( INFO.NE.0 )
513: $ GO TO 70
514: *
515: * Multiply B by reciprocals of singular values
516: *
517: THR = MAX( RCOND*S( 1 ), SFMIN )
518: IF( RCOND.LT.ZERO )
519: $ THR = MAX( EPS*S( 1 ), SFMIN )
520: RANK = 0
521: DO 10 I = 1, N
522: IF( S( I ).GT.THR ) THEN
523: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
524: RANK = RANK + 1
525: ELSE
526: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
527: END IF
528: 10 CONTINUE
529: *
530: * Multiply B by right singular vectors
531: * (CWorkspace: need N, prefer N*NRHS)
532: * (RWorkspace: none)
533: *
534: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
535: CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
536: $ CZERO, WORK, LDB )
537: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
538: ELSE IF( NRHS.GT.1 ) THEN
539: CHUNK = LWORK / N
540: DO 20 I = 1, NRHS, CHUNK
541: BL = MIN( NRHS-I+1, CHUNK )
542: CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
543: $ LDB, CZERO, WORK, N )
544: CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
545: 20 CONTINUE
546: ELSE
547: CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
548: CALL ZCOPY( N, WORK, 1, B, 1 )
549: END IF
550: *
551: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
552: $ THEN
553: *
554: * Underdetermined case, M much less than N
555: *
556: * Path 2a - underdetermined, with many more columns than rows
557: * and sufficient workspace for an efficient algorithm
558: *
559: LDWORK = M
560: IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
561: $ LDWORK = LDA
562: ITAU = 1
563: IWORK = M + 1
564: *
565: * Compute A=L*Q
566: * (CWorkspace: need 2*M, prefer M+M*NB)
567: * (RWorkspace: none)
568: *
569: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
570: $ LWORK-IWORK+1, INFO )
571: IL = IWORK
572: *
573: * Copy L to WORK(IL), zeroing out above it
574: *
575: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
576: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
577: $ LDWORK )
578: IE = 1
579: ITAUQ = IL + LDWORK*M
580: ITAUP = ITAUQ + M
581: IWORK = ITAUP + M
582: *
583: * Bidiagonalize L in WORK(IL)
584: * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
585: * (RWorkspace: need M)
586: *
587: CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
588: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
589: $ LWORK-IWORK+1, INFO )
590: *
591: * Multiply B by transpose of left bidiagonalizing vectors of L
592: * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
593: * (RWorkspace: none)
594: *
595: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
596: $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
597: $ LWORK-IWORK+1, INFO )
598: *
599: * Generate right bidiagonalizing vectors of R in WORK(IL)
600: * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
601: * (RWorkspace: none)
602: *
603: CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
604: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
605: IRWORK = IE + M
606: *
607: * Perform bidiagonal QR iteration, computing right singular
608: * vectors of L in WORK(IL) and multiplying B by transpose of
609: * left singular vectors
610: * (CWorkspace: need M*M)
611: * (RWorkspace: need BDSPAC)
612: *
613: CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
614: $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
615: IF( INFO.NE.0 )
616: $ GO TO 70
617: *
618: * Multiply B by reciprocals of singular values
619: *
620: THR = MAX( RCOND*S( 1 ), SFMIN )
621: IF( RCOND.LT.ZERO )
622: $ THR = MAX( EPS*S( 1 ), SFMIN )
623: RANK = 0
624: DO 30 I = 1, M
625: IF( S( I ).GT.THR ) THEN
626: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
627: RANK = RANK + 1
628: ELSE
629: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
630: END IF
631: 30 CONTINUE
632: IWORK = IL + M*LDWORK
633: *
634: * Multiply B by right singular vectors of L in WORK(IL)
635: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
636: * (RWorkspace: none)
637: *
638: IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
639: CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
640: $ B, LDB, CZERO, WORK( IWORK ), LDB )
641: CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
642: ELSE IF( NRHS.GT.1 ) THEN
643: CHUNK = ( LWORK-IWORK+1 ) / M
644: DO 40 I = 1, NRHS, CHUNK
645: BL = MIN( NRHS-I+1, CHUNK )
646: CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
647: $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
648: CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
649: $ LDB )
650: 40 CONTINUE
651: ELSE
652: CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
653: $ 1, CZERO, WORK( IWORK ), 1 )
654: CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
655: END IF
656: *
657: * Zero out below first M rows of B
658: *
659: CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
660: IWORK = ITAU + M
661: *
662: * Multiply transpose(Q) by B
663: * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
664: * (RWorkspace: none)
665: *
666: CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
667: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
668: *
669: ELSE
670: *
671: * Path 2 - remaining underdetermined cases
672: *
673: IE = 1
674: ITAUQ = 1
675: ITAUP = ITAUQ + M
676: IWORK = ITAUP + M
677: *
678: * Bidiagonalize A
679: * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
680: * (RWorkspace: need N)
681: *
682: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
683: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
684: $ INFO )
685: *
686: * Multiply B by transpose of left bidiagonalizing vectors
687: * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
688: * (RWorkspace: none)
689: *
690: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
691: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
692: *
693: * Generate right bidiagonalizing vectors in A
694: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
695: * (RWorkspace: none)
696: *
697: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
698: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
699: IRWORK = IE + M
700: *
701: * Perform bidiagonal QR iteration,
702: * computing right singular vectors of A in A and
703: * multiplying B by transpose of left singular vectors
704: * (CWorkspace: none)
705: * (RWorkspace: need BDSPAC)
706: *
707: CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
708: $ 1, B, LDB, RWORK( IRWORK ), INFO )
709: IF( INFO.NE.0 )
710: $ GO TO 70
711: *
712: * Multiply B by reciprocals of singular values
713: *
714: THR = MAX( RCOND*S( 1 ), SFMIN )
715: IF( RCOND.LT.ZERO )
716: $ THR = MAX( EPS*S( 1 ), SFMIN )
717: RANK = 0
718: DO 50 I = 1, M
719: IF( S( I ).GT.THR ) THEN
720: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
721: RANK = RANK + 1
722: ELSE
723: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
724: END IF
725: 50 CONTINUE
726: *
727: * Multiply B by right singular vectors of A
728: * (CWorkspace: need N, prefer N*NRHS)
729: * (RWorkspace: none)
730: *
731: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
732: CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
733: $ CZERO, WORK, LDB )
734: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
735: ELSE IF( NRHS.GT.1 ) THEN
736: CHUNK = LWORK / N
737: DO 60 I = 1, NRHS, CHUNK
738: BL = MIN( NRHS-I+1, CHUNK )
739: CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
740: $ LDB, CZERO, WORK, N )
741: CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
742: 60 CONTINUE
743: ELSE
744: CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
745: CALL ZCOPY( N, WORK, 1, B, 1 )
746: END IF
747: END IF
748: *
749: * Undo scaling
750: *
751: IF( IASCL.EQ.1 ) THEN
752: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
753: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
754: $ INFO )
755: ELSE IF( IASCL.EQ.2 ) THEN
756: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
757: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
758: $ INFO )
759: END IF
760: IF( IBSCL.EQ.1 ) THEN
761: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
762: ELSE IF( IBSCL.EQ.2 ) THEN
763: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
764: END IF
765: 70 CONTINUE
766: WORK( 1 ) = MAXWRK
767: RETURN
768: *
769: * End of ZGELSS
770: *
771: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelss.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack