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: *> \date June 2017
214: *
215: *> \ingroup complex16GEsolve
216: *
217: *> \par Contributors:
218: * ==================
219: *>
220: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
221: *> California at Berkeley, USA \n
222: *> Osni Marques, LBNL/NERSC, USA \n
223: *
224: * =====================================================================
225: SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
226: $ WORK, LWORK, RWORK, IWORK, INFO )
227: *
228: * -- LAPACK driver routine (version 3.7.1) --
229: * -- LAPACK is a software package provided by Univ. of Tennessee, --
230: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231: * June 2017
232: *
233: * .. Scalar Arguments ..
234: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
235: DOUBLE PRECISION RCOND
236: * ..
237: * .. Array Arguments ..
238: INTEGER IWORK( * )
239: DOUBLE PRECISION RWORK( * ), S( * )
240: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
241: * ..
242: *
243: * =====================================================================
244: *
245: * .. Parameters ..
246: DOUBLE PRECISION ZERO, ONE, TWO
247: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
248: COMPLEX*16 CZERO
249: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
250: * ..
251: * .. Local Scalars ..
252: LOGICAL LQUERY
253: INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
254: $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
255: $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
256: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
257: * ..
258: * .. External Subroutines ..
259: EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF,
260: $ ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR,
261: $ ZUNMLQ, ZUNMQR
262: * ..
263: * .. External Functions ..
264: INTEGER ILAENV
265: DOUBLE PRECISION DLAMCH, ZLANGE
266: EXTERNAL ILAENV, DLAMCH, ZLANGE
267: * ..
268: * .. Intrinsic Functions ..
269: INTRINSIC INT, LOG, MAX, MIN, DBLE
270: * ..
271: * .. Executable Statements ..
272: *
273: * Test the input arguments.
274: *
275: INFO = 0
276: MINMN = MIN( M, N )
277: MAXMN = MAX( M, N )
278: LQUERY = ( LWORK.EQ.-1 )
279: IF( M.LT.0 ) THEN
280: INFO = -1
281: ELSE IF( N.LT.0 ) THEN
282: INFO = -2
283: ELSE IF( NRHS.LT.0 ) THEN
284: INFO = -3
285: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
286: INFO = -5
287: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
288: INFO = -7
289: END IF
290: *
291: * Compute workspace.
292: * (Note: Comments in the code beginning "Workspace:" describe the
293: * minimal amount of workspace needed at that point in the code,
294: * as well as the preferred amount for good performance.
295: * NB refers to the optimal block size for the immediately
296: * following subroutine, as returned by ILAENV.)
297: *
298: IF( INFO.EQ.0 ) THEN
299: MINWRK = 1
300: MAXWRK = 1
301: LIWORK = 1
302: LRWORK = 1
303: IF( MINMN.GT.0 ) THEN
304: SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 )
305: MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 )
306: NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) /
307: $ LOG( TWO ) ) + 1, 0 )
308: LIWORK = 3*MINMN*NLVL + 11*MINMN
309: MM = M
310: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
311: *
312: * Path 1a - overdetermined, with many more rows than
313: * columns.
314: *
315: MM = N
316: MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N,
317: $ -1, -1 ) )
318: MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M,
319: $ NRHS, N, -1 ) )
320: END IF
321: IF( M.GE.N ) THEN
322: *
323: * Path 1 - overdetermined or exactly determined.
324: *
325: LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
326: $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
327: MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
328: $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
329: MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
330: $ 'QLC', MM, NRHS, N, -1 ) )
331: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
332: $ 'ZUNMBR', 'PLN', N, NRHS, N, -1 ) )
333: MAXWRK = MAX( MAXWRK, 2*N + N*NRHS )
334: MINWRK = MAX( 2*N + MM, 2*N + N*NRHS )
335: END IF
336: IF( N.GT.M ) THEN
337: LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS +
338: $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS )
339: IF( N.GE.MNTHR ) THEN
340: *
341: * Path 2a - underdetermined, with many more columns
342: * than rows.
343: *
344: MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
345: $ -1 )
346: MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1,
347: $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
348: MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1,
349: $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
350: MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1,
351: $ 'ZUNMLQ', 'LC', N, NRHS, M, -1 ) )
352: IF( NRHS.GT.1 ) THEN
353: MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
354: ELSE
355: MAXWRK = MAX( MAXWRK, M*M + 2*M )
356: END IF
357: MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS )
358: ! XXX: Ensure the Path 2a case below is triggered. The workspace
359: ! calculation should use queries for all routines eventually.
360: MAXWRK = MAX( MAXWRK,
361: $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
362: ELSE
363: *
364: * Path 2 - underdetermined.
365: *
366: MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
367: $ N, -1, -1 )
368: MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
369: $ 'QLC', M, NRHS, M, -1 ) )
370: MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR',
371: $ 'PLN', N, NRHS, M, -1 ) )
372: MAXWRK = MAX( MAXWRK, 2*M + M*NRHS )
373: END IF
374: MINWRK = MAX( 2*M + N, 2*M + M*NRHS )
375: END IF
376: END IF
377: MINWRK = MIN( MINWRK, MAXWRK )
378: WORK( 1 ) = MAXWRK
379: IWORK( 1 ) = LIWORK
380: RWORK( 1 ) = LRWORK
381: *
382: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
383: INFO = -12
384: END IF
385: END IF
386: *
387: IF( INFO.NE.0 ) THEN
388: CALL XERBLA( 'ZGELSD', -INFO )
389: RETURN
390: ELSE IF( LQUERY ) THEN
391: RETURN
392: END IF
393: *
394: * Quick return if possible.
395: *
396: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
397: RANK = 0
398: RETURN
399: END IF
400: *
401: * Get machine parameters.
402: *
403: EPS = DLAMCH( 'P' )
404: SFMIN = DLAMCH( 'S' )
405: SMLNUM = SFMIN / EPS
406: BIGNUM = ONE / SMLNUM
407: CALL DLABAD( SMLNUM, BIGNUM )
408: *
409: * Scale A if max entry outside range [SMLNUM,BIGNUM].
410: *
411: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
412: IASCL = 0
413: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
414: *
415: * Scale matrix norm up to SMLNUM
416: *
417: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
418: IASCL = 1
419: ELSE IF( ANRM.GT.BIGNUM ) THEN
420: *
421: * Scale matrix norm down to BIGNUM.
422: *
423: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
424: IASCL = 2
425: ELSE IF( ANRM.EQ.ZERO ) THEN
426: *
427: * Matrix all zero. Return zero solution.
428: *
429: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
430: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
431: RANK = 0
432: GO TO 10
433: END IF
434: *
435: * Scale B if max entry outside range [SMLNUM,BIGNUM].
436: *
437: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
438: IBSCL = 0
439: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
440: *
441: * Scale matrix norm up to SMLNUM.
442: *
443: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
444: IBSCL = 1
445: ELSE IF( BNRM.GT.BIGNUM ) THEN
446: *
447: * Scale matrix norm down to BIGNUM.
448: *
449: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
450: IBSCL = 2
451: END IF
452: *
453: * If M < N make sure B(M+1:N,:) = 0
454: *
455: IF( M.LT.N )
456: $ CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
457: *
458: * Overdetermined case.
459: *
460: IF( M.GE.N ) THEN
461: *
462: * Path 1 - overdetermined or exactly determined.
463: *
464: MM = M
465: IF( M.GE.MNTHR ) THEN
466: *
467: * Path 1a - overdetermined, with many more rows than columns
468: *
469: MM = N
470: ITAU = 1
471: NWORK = ITAU + N
472: *
473: * Compute A=Q*R.
474: * (RWorkspace: need N)
475: * (CWorkspace: need N, prefer N*NB)
476: *
477: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
478: $ LWORK-NWORK+1, INFO )
479: *
480: * Multiply B by transpose(Q).
481: * (RWorkspace: need N)
482: * (CWorkspace: need NRHS, prefer NRHS*NB)
483: *
484: CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
485: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
486: *
487: * Zero out below R.
488: *
489: IF( N.GT.1 ) THEN
490: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
491: $ LDA )
492: END IF
493: END IF
494: *
495: ITAUQ = 1
496: ITAUP = ITAUQ + N
497: NWORK = ITAUP + N
498: IE = 1
499: NRWORK = IE + N
500: *
501: * Bidiagonalize R in A.
502: * (RWorkspace: need N)
503: * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
504: *
505: CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
506: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
507: $ INFO )
508: *
509: * Multiply B by transpose of left bidiagonalizing vectors of R.
510: * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
511: *
512: CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
513: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
514: *
515: * Solve the bidiagonal least squares problem.
516: *
517: CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
518: $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
519: $ IWORK, INFO )
520: IF( INFO.NE.0 ) THEN
521: GO TO 10
522: END IF
523: *
524: * Multiply B by right bidiagonalizing vectors of R.
525: *
526: CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
527: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
528: *
529: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
530: $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
531: *
532: * Path 2a - underdetermined, with many more columns than rows
533: * and sufficient workspace for an efficient algorithm.
534: *
535: LDWORK = M
536: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
537: $ M*LDA+M+M*NRHS ) )LDWORK = LDA
538: ITAU = 1
539: NWORK = M + 1
540: *
541: * Compute A=L*Q.
542: * (CWorkspace: need 2*M, prefer M+M*NB)
543: *
544: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
545: $ LWORK-NWORK+1, INFO )
546: IL = NWORK
547: *
548: * Copy L to WORK(IL), zeroing out above its diagonal.
549: *
550: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
551: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
552: $ LDWORK )
553: ITAUQ = IL + LDWORK*M
554: ITAUP = ITAUQ + M
555: NWORK = ITAUP + M
556: IE = 1
557: NRWORK = IE + M
558: *
559: * Bidiagonalize L in WORK(IL).
560: * (RWorkspace: need M)
561: * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
562: *
563: CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
564: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
565: $ LWORK-NWORK+1, INFO )
566: *
567: * Multiply B by transpose of left bidiagonalizing vectors of L.
568: * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
569: *
570: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
571: $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
572: $ LWORK-NWORK+1, INFO )
573: *
574: * Solve the bidiagonal least squares problem.
575: *
576: CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
577: $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
578: $ IWORK, INFO )
579: IF( INFO.NE.0 ) THEN
580: GO TO 10
581: END IF
582: *
583: * Multiply B by right bidiagonalizing vectors of L.
584: *
585: CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
586: $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
587: $ LWORK-NWORK+1, INFO )
588: *
589: * Zero out below first M rows of B.
590: *
591: CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
592: NWORK = ITAU + M
593: *
594: * Multiply transpose(Q) by B.
595: * (CWorkspace: need NRHS, prefer NRHS*NB)
596: *
597: CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
598: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
599: *
600: ELSE
601: *
602: * Path 2 - remaining underdetermined cases.
603: *
604: ITAUQ = 1
605: ITAUP = ITAUQ + M
606: NWORK = ITAUP + M
607: IE = 1
608: NRWORK = IE + M
609: *
610: * Bidiagonalize A.
611: * (RWorkspace: need M)
612: * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
613: *
614: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
615: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
616: $ INFO )
617: *
618: * Multiply B by transpose of left bidiagonalizing vectors.
619: * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
620: *
621: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
622: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
623: *
624: * Solve the bidiagonal least squares problem.
625: *
626: CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
627: $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
628: $ IWORK, INFO )
629: IF( INFO.NE.0 ) THEN
630: GO TO 10
631: END IF
632: *
633: * Multiply B by right bidiagonalizing vectors of A.
634: *
635: CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
636: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
637: *
638: END IF
639: *
640: * Undo scaling.
641: *
642: IF( IASCL.EQ.1 ) THEN
643: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
644: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
645: $ INFO )
646: ELSE IF( IASCL.EQ.2 ) THEN
647: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
648: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
649: $ INFO )
650: END IF
651: IF( IBSCL.EQ.1 ) THEN
652: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
653: ELSE IF( IBSCL.EQ.2 ) THEN
654: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
655: END IF
656: *
657: 10 CONTINUE
658: WORK( 1 ) = MAXWRK
659: IWORK( 1 ) = LIWORK
660: RWORK( 1 ) = LRWORK
661: RETURN
662: *
663: * End of ZGELSD
664: *
665: 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