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