1: SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
2: $ WORK, LWORK, RWORK, INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
11: DOUBLE PRECISION RCOND
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION RWORK( * ), S( * )
15: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZGELSS computes the minimum norm solution to a complex linear
22: * least squares problem:
23: *
24: * Minimize 2-norm(| b - A*x |).
25: *
26: * using the singular value decomposition (SVD) of A. A is an M-by-N
27: * matrix which may be rank-deficient.
28: *
29: * Several right hand side vectors b and solution vectors x can be
30: * handled in a single call; they are stored as the columns of the
31: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
32: * X.
33: *
34: * The effective rank of A is determined by treating as zero those
35: * singular values which are less than RCOND times the largest singular
36: * value.
37: *
38: * Arguments
39: * =========
40: *
41: * M (input) INTEGER
42: * The number of rows of the matrix A. M >= 0.
43: *
44: * N (input) INTEGER
45: * The number of columns of the matrix A. N >= 0.
46: *
47: * NRHS (input) INTEGER
48: * The number of right hand sides, i.e., the number of columns
49: * of the matrices B and X. NRHS >= 0.
50: *
51: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
52: * On entry, the M-by-N matrix A.
53: * On exit, the first min(m,n) rows of A are overwritten with
54: * its right singular vectors, stored rowwise.
55: *
56: * LDA (input) INTEGER
57: * The leading dimension of the array A. LDA >= max(1,M).
58: *
59: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
60: * On entry, the M-by-NRHS right hand side matrix B.
61: * On exit, B is overwritten by the N-by-NRHS solution matrix X.
62: * If m >= n and RANK = n, the residual sum-of-squares for
63: * the solution in the i-th column is given by the sum of
64: * squares of the modulus of elements n+1:m in that column.
65: *
66: * LDB (input) INTEGER
67: * The leading dimension of the array B. LDB >= max(1,M,N).
68: *
69: * S (output) DOUBLE PRECISION array, dimension (min(M,N))
70: * The singular values of A in decreasing order.
71: * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
72: *
73: * RCOND (input) DOUBLE PRECISION
74: * RCOND is used to determine the effective rank of A.
75: * Singular values S(i) <= RCOND*S(1) are treated as zero.
76: * If RCOND < 0, machine precision is used instead.
77: *
78: * RANK (output) INTEGER
79: * The effective rank of A, i.e., the number of singular values
80: * which are greater than RCOND*S(1).
81: *
82: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
83: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84: *
85: * LWORK (input) INTEGER
86: * The dimension of the array WORK. LWORK >= 1, and also:
87: * LWORK >= 2*min(M,N) + max(M,N,NRHS)
88: * For good performance, LWORK should generally be larger.
89: *
90: * If LWORK = -1, then a workspace query is assumed; the routine
91: * only calculates the optimal size of the WORK array, returns
92: * this value as the first entry of the WORK array, and no error
93: * message related to LWORK is issued by XERBLA.
94: *
95: * RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
96: *
97: * INFO (output) INTEGER
98: * = 0: successful exit
99: * < 0: if INFO = -i, the i-th argument had an illegal value.
100: * > 0: the algorithm for computing the SVD failed to converge;
101: * if INFO = i, i off-diagonal elements of an intermediate
102: * bidiagonal form did not converge to zero.
103: *
104: * =====================================================================
105: *
106: * .. Parameters ..
107: DOUBLE PRECISION ZERO, ONE
108: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
109: COMPLEX*16 CZERO, CONE
110: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
111: $ CONE = ( 1.0D+0, 0.0D+0 ) )
112: * ..
113: * .. Local Scalars ..
114: LOGICAL LQUERY
115: INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
116: $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
117: $ MAXWRK, MINMN, MINWRK, MM, MNTHR
118: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
119: * ..
120: * .. Local Arrays ..
121: COMPLEX*16 VDUM( 1 )
122: * ..
123: * .. External Subroutines ..
124: EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY,
125: $ ZDRSCL, ZGEBRD, ZGELQF, ZGEMM, ZGEMV, ZGEQRF,
126: $ ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNMBR, ZUNMLQ,
127: $ ZUNMQR
128: * ..
129: * .. External Functions ..
130: INTEGER ILAENV
131: DOUBLE PRECISION DLAMCH, ZLANGE
132: EXTERNAL ILAENV, DLAMCH, ZLANGE
133: * ..
134: * .. Intrinsic Functions ..
135: INTRINSIC MAX, MIN
136: * ..
137: * .. Executable Statements ..
138: *
139: * Test the input arguments
140: *
141: INFO = 0
142: MINMN = MIN( M, N )
143: MAXMN = MAX( M, N )
144: LQUERY = ( LWORK.EQ.-1 )
145: IF( M.LT.0 ) THEN
146: INFO = -1
147: ELSE IF( N.LT.0 ) THEN
148: INFO = -2
149: ELSE IF( NRHS.LT.0 ) THEN
150: INFO = -3
151: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
152: INFO = -5
153: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
154: INFO = -7
155: END IF
156: *
157: * Compute workspace
158: * (Note: Comments in the code beginning "Workspace:" describe the
159: * minimal amount of workspace needed at that point in the code,
160: * as well as the preferred amount for good performance.
161: * CWorkspace refers to complex workspace, and RWorkspace refers
162: * to real workspace. NB refers to the optimal block size for the
163: * immediately following subroutine, as returned by ILAENV.)
164: *
165: IF( INFO.EQ.0 ) THEN
166: MINWRK = 1
167: MAXWRK = 1
168: IF( MINMN.GT.0 ) THEN
169: MM = M
170: MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
171: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
172: *
173: * Path 1a - overdetermined, with many more rows than
174: * columns
175: *
176: MM = N
177: MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
178: $ N, -1, -1 ) )
179: MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
180: $ M, NRHS, N, -1 ) )
181: END IF
182: IF( M.GE.N ) THEN
183: *
184: * Path 1 - overdetermined or exactly determined
185: *
186: MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
187: $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
188: MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
189: $ 'QLC', MM, NRHS, N, -1 ) )
190: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
191: $ 'ZUNGBR', 'P', N, N, N, -1 ) )
192: MAXWRK = MAX( MAXWRK, N*NRHS )
193: MINWRK = 2*N + MAX( NRHS, M )
194: END IF
195: IF( N.GT.M ) THEN
196: MINWRK = 2*M + MAX( NRHS, N )
197: IF( N.GE.MNTHR ) THEN
198: *
199: * Path 2a - underdetermined, with many more columns
200: * than rows
201: *
202: MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
203: $ -1 )
204: MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
205: $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
206: MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
207: $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
208: MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
209: $ 'ZUNGBR', 'P', M, M, M, -1 ) )
210: IF( NRHS.GT.1 ) THEN
211: MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
212: ELSE
213: MAXWRK = MAX( MAXWRK, M*M + 2*M )
214: END IF
215: MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ',
216: $ 'LC', N, NRHS, M, -1 ) )
217: ELSE
218: *
219: * Path 2 - underdetermined
220: *
221: MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
222: $ N, -1, -1 )
223: MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
224: $ 'QLC', M, NRHS, M, -1 ) )
225: MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR',
226: $ 'P', M, N, M, -1 ) )
227: MAXWRK = MAX( MAXWRK, N*NRHS )
228: END IF
229: END IF
230: MAXWRK = MAX( MINWRK, MAXWRK )
231: END IF
232: WORK( 1 ) = MAXWRK
233: *
234: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
235: $ INFO = -12
236: END IF
237: *
238: IF( INFO.NE.0 ) THEN
239: CALL XERBLA( 'ZGELSS', -INFO )
240: RETURN
241: ELSE IF( LQUERY ) THEN
242: RETURN
243: END IF
244: *
245: * Quick return if possible
246: *
247: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
248: RANK = 0
249: RETURN
250: END IF
251: *
252: * Get machine parameters
253: *
254: EPS = DLAMCH( 'P' )
255: SFMIN = DLAMCH( 'S' )
256: SMLNUM = SFMIN / EPS
257: BIGNUM = ONE / SMLNUM
258: CALL DLABAD( SMLNUM, BIGNUM )
259: *
260: * Scale A if max element outside range [SMLNUM,BIGNUM]
261: *
262: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
263: IASCL = 0
264: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
265: *
266: * Scale matrix norm up to SMLNUM
267: *
268: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
269: IASCL = 1
270: ELSE IF( ANRM.GT.BIGNUM ) THEN
271: *
272: * Scale matrix norm down to BIGNUM
273: *
274: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
275: IASCL = 2
276: ELSE IF( ANRM.EQ.ZERO ) THEN
277: *
278: * Matrix all zero. Return zero solution.
279: *
280: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
281: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
282: RANK = 0
283: GO TO 70
284: END IF
285: *
286: * Scale B if max element outside range [SMLNUM,BIGNUM]
287: *
288: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
289: IBSCL = 0
290: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
291: *
292: * Scale matrix norm up to SMLNUM
293: *
294: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
295: IBSCL = 1
296: ELSE IF( BNRM.GT.BIGNUM ) THEN
297: *
298: * Scale matrix norm down to BIGNUM
299: *
300: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
301: IBSCL = 2
302: END IF
303: *
304: * Overdetermined case
305: *
306: IF( M.GE.N ) THEN
307: *
308: * Path 1 - overdetermined or exactly determined
309: *
310: MM = M
311: IF( M.GE.MNTHR ) THEN
312: *
313: * Path 1a - overdetermined, with many more rows than columns
314: *
315: MM = N
316: ITAU = 1
317: IWORK = ITAU + N
318: *
319: * Compute A=Q*R
320: * (CWorkspace: need 2*N, prefer N+N*NB)
321: * (RWorkspace: none)
322: *
323: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
324: $ LWORK-IWORK+1, INFO )
325: *
326: * Multiply B by transpose(Q)
327: * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
328: * (RWorkspace: none)
329: *
330: CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
331: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
332: *
333: * Zero out below R
334: *
335: IF( N.GT.1 )
336: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
337: $ LDA )
338: END IF
339: *
340: IE = 1
341: ITAUQ = 1
342: ITAUP = ITAUQ + N
343: IWORK = ITAUP + N
344: *
345: * Bidiagonalize R in A
346: * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
347: * (RWorkspace: need N)
348: *
349: CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
350: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
351: $ INFO )
352: *
353: * Multiply B by transpose of left bidiagonalizing vectors of R
354: * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
355: * (RWorkspace: none)
356: *
357: CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
358: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
359: *
360: * Generate right bidiagonalizing vectors of R in A
361: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
362: * (RWorkspace: none)
363: *
364: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
365: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
366: IRWORK = IE + N
367: *
368: * Perform bidiagonal QR iteration
369: * multiply B by transpose of left singular vectors
370: * compute right singular vectors in A
371: * (CWorkspace: none)
372: * (RWorkspace: need BDSPAC)
373: *
374: CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
375: $ 1, B, LDB, RWORK( IRWORK ), INFO )
376: IF( INFO.NE.0 )
377: $ GO TO 70
378: *
379: * Multiply B by reciprocals of singular values
380: *
381: THR = MAX( RCOND*S( 1 ), SFMIN )
382: IF( RCOND.LT.ZERO )
383: $ THR = MAX( EPS*S( 1 ), SFMIN )
384: RANK = 0
385: DO 10 I = 1, N
386: IF( S( I ).GT.THR ) THEN
387: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
388: RANK = RANK + 1
389: ELSE
390: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
391: END IF
392: 10 CONTINUE
393: *
394: * Multiply B by right singular vectors
395: * (CWorkspace: need N, prefer N*NRHS)
396: * (RWorkspace: none)
397: *
398: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
399: CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
400: $ CZERO, WORK, LDB )
401: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
402: ELSE IF( NRHS.GT.1 ) THEN
403: CHUNK = LWORK / N
404: DO 20 I = 1, NRHS, CHUNK
405: BL = MIN( NRHS-I+1, CHUNK )
406: CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
407: $ LDB, CZERO, WORK, N )
408: CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
409: 20 CONTINUE
410: ELSE
411: CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
412: CALL ZCOPY( N, WORK, 1, B, 1 )
413: END IF
414: *
415: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
416: $ THEN
417: *
418: * Underdetermined case, M much less than N
419: *
420: * Path 2a - underdetermined, with many more columns than rows
421: * and sufficient workspace for an efficient algorithm
422: *
423: LDWORK = M
424: IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
425: $ LDWORK = LDA
426: ITAU = 1
427: IWORK = M + 1
428: *
429: * Compute A=L*Q
430: * (CWorkspace: need 2*M, prefer M+M*NB)
431: * (RWorkspace: none)
432: *
433: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
434: $ LWORK-IWORK+1, INFO )
435: IL = IWORK
436: *
437: * Copy L to WORK(IL), zeroing out above it
438: *
439: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
440: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
441: $ LDWORK )
442: IE = 1
443: ITAUQ = IL + LDWORK*M
444: ITAUP = ITAUQ + M
445: IWORK = ITAUP + M
446: *
447: * Bidiagonalize L in WORK(IL)
448: * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
449: * (RWorkspace: need M)
450: *
451: CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
452: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
453: $ LWORK-IWORK+1, INFO )
454: *
455: * Multiply B by transpose of left bidiagonalizing vectors of L
456: * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
457: * (RWorkspace: none)
458: *
459: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
460: $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
461: $ LWORK-IWORK+1, INFO )
462: *
463: * Generate right bidiagonalizing vectors of R in WORK(IL)
464: * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
465: * (RWorkspace: none)
466: *
467: CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
468: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
469: IRWORK = IE + M
470: *
471: * Perform bidiagonal QR iteration, computing right singular
472: * vectors of L in WORK(IL) and multiplying B by transpose of
473: * left singular vectors
474: * (CWorkspace: need M*M)
475: * (RWorkspace: need BDSPAC)
476: *
477: CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
478: $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
479: IF( INFO.NE.0 )
480: $ GO TO 70
481: *
482: * Multiply B by reciprocals of singular values
483: *
484: THR = MAX( RCOND*S( 1 ), SFMIN )
485: IF( RCOND.LT.ZERO )
486: $ THR = MAX( EPS*S( 1 ), SFMIN )
487: RANK = 0
488: DO 30 I = 1, M
489: IF( S( I ).GT.THR ) THEN
490: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
491: RANK = RANK + 1
492: ELSE
493: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
494: END IF
495: 30 CONTINUE
496: IWORK = IL + M*LDWORK
497: *
498: * Multiply B by right singular vectors of L in WORK(IL)
499: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
500: * (RWorkspace: none)
501: *
502: IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
503: CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
504: $ B, LDB, CZERO, WORK( IWORK ), LDB )
505: CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
506: ELSE IF( NRHS.GT.1 ) THEN
507: CHUNK = ( LWORK-IWORK+1 ) / M
508: DO 40 I = 1, NRHS, CHUNK
509: BL = MIN( NRHS-I+1, CHUNK )
510: CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
511: $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
512: CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
513: $ LDB )
514: 40 CONTINUE
515: ELSE
516: CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
517: $ 1, CZERO, WORK( IWORK ), 1 )
518: CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
519: END IF
520: *
521: * Zero out below first M rows of B
522: *
523: CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
524: IWORK = ITAU + M
525: *
526: * Multiply transpose(Q) by B
527: * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
528: * (RWorkspace: none)
529: *
530: CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
531: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
532: *
533: ELSE
534: *
535: * Path 2 - remaining underdetermined cases
536: *
537: IE = 1
538: ITAUQ = 1
539: ITAUP = ITAUQ + M
540: IWORK = ITAUP + M
541: *
542: * Bidiagonalize A
543: * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
544: * (RWorkspace: need N)
545: *
546: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
547: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
548: $ INFO )
549: *
550: * Multiply B by transpose of left bidiagonalizing vectors
551: * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
552: * (RWorkspace: none)
553: *
554: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
555: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
556: *
557: * Generate right bidiagonalizing vectors in A
558: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
559: * (RWorkspace: none)
560: *
561: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
562: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
563: IRWORK = IE + M
564: *
565: * Perform bidiagonal QR iteration,
566: * computing right singular vectors of A in A and
567: * multiplying B by transpose of left singular vectors
568: * (CWorkspace: none)
569: * (RWorkspace: need BDSPAC)
570: *
571: CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
572: $ 1, B, LDB, RWORK( IRWORK ), INFO )
573: IF( INFO.NE.0 )
574: $ GO TO 70
575: *
576: * Multiply B by reciprocals of singular values
577: *
578: THR = MAX( RCOND*S( 1 ), SFMIN )
579: IF( RCOND.LT.ZERO )
580: $ THR = MAX( EPS*S( 1 ), SFMIN )
581: RANK = 0
582: DO 50 I = 1, M
583: IF( S( I ).GT.THR ) THEN
584: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
585: RANK = RANK + 1
586: ELSE
587: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
588: END IF
589: 50 CONTINUE
590: *
591: * Multiply B by right singular vectors of A
592: * (CWorkspace: need N, prefer N*NRHS)
593: * (RWorkspace: none)
594: *
595: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
596: CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
597: $ CZERO, WORK, LDB )
598: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
599: ELSE IF( NRHS.GT.1 ) THEN
600: CHUNK = LWORK / N
601: DO 60 I = 1, NRHS, CHUNK
602: BL = MIN( NRHS-I+1, CHUNK )
603: CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
604: $ LDB, CZERO, WORK, N )
605: CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
606: 60 CONTINUE
607: ELSE
608: CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
609: CALL ZCOPY( N, WORK, 1, B, 1 )
610: END IF
611: END IF
612: *
613: * Undo scaling
614: *
615: IF( IASCL.EQ.1 ) THEN
616: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
617: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
618: $ INFO )
619: ELSE IF( IASCL.EQ.2 ) THEN
620: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
621: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
622: $ INFO )
623: END IF
624: IF( IBSCL.EQ.1 ) THEN
625: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
626: ELSE IF( IBSCL.EQ.2 ) THEN
627: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
628: END IF
629: 70 CONTINUE
630: WORK( 1 ) = MAXWRK
631: RETURN
632: *
633: * End of ZGELSS
634: *
635: 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