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