OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dgesdd.f
Go to the documentation of this file.
1*> \brief \b DGESDD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGESDD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22* WORK, LWORK, IWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBZ
26* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
27* ..
28* .. Array Arguments ..
29* INTEGER IWORK( * )
30* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
31* $ VT( LDVT, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DGESDD computes the singular value decomposition (SVD) of a real
41*> M-by-N matrix A, optionally computing the left and right singular
42*> vectors. If singular vectors are desired, it uses a
43*> divide-and-conquer algorithm.
44*>
45*> The SVD is written
46*>
47*> A = U * SIGMA * transpose(V)
48*>
49*> where SIGMA is an M-by-N matrix which is zero except for its
50*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
51*> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
52*> are the singular values of A; they are real and non-negative, and
53*> are returned in descending order. The first min(m,n) columns of
54*> U and V are the left and right singular vectors of A.
55*>
56*> Note that the routine returns VT = V**T, not V.
57*>
58*> The divide and conquer algorithm makes very mild assumptions about
59*> floating point arithmetic. It will work on machines with a guard
60*> digit in add/subtract, or on those binary machines without guard
61*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
62*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
63*> without guard digits, but we know of none.
64*> \endverbatim
65*
66* Arguments:
67* ==========
68*
69*> \param[in] JOBZ
70*> \verbatim
71*> JOBZ is CHARACTER*1
72*> Specifies options for computing all or part of the matrix U:
73*> = 'A': all M columns of U and all N rows of V**T are
74*> returned in the arrays U and VT;
75*> = 'S': the first min(M,N) columns of U and the first
76*> min(M,N) rows of V**T are returned in the arrays U
77*> and VT;
78*> = 'O': If M >= N, the first N columns of U are overwritten
79*> on the array A and all rows of V**T are returned in
80*> the array VT;
81*> otherwise, all columns of U are returned in the
82*> array U and the first M rows of V**T are overwritten
83*> in the array A;
84*> = 'N': no columns of U or rows of V**T are computed.
85*> \endverbatim
86*>
87*> \param[in] M
88*> \verbatim
89*> M is INTEGER
90*> The number of rows of the input matrix A. M >= 0.
91*> \endverbatim
92*>
93*> \param[in] N
94*> \verbatim
95*> N is INTEGER
96*> The number of columns of the input matrix A. N >= 0.
97*> \endverbatim
98*>
99*> \param[in,out] A
100*> \verbatim
101*> A is DOUBLE PRECISION array, dimension (LDA,N)
102*> On entry, the M-by-N matrix A.
103*> On exit,
104*> if JOBZ = 'O', A is overwritten with the first N columns
105*> of U (the left singular vectors, stored
106*> columnwise) if M >= N;
107*> A is overwritten with the first M rows
108*> of V**T (the right singular vectors, stored
109*> rowwise) otherwise.
110*> if JOBZ .ne. 'O', the contents of A are destroyed.
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*> LDA is INTEGER
116*> The leading dimension of the array A. LDA >= max(1,M).
117*> \endverbatim
118*>
119*> \param[out] S
120*> \verbatim
121*> S is DOUBLE PRECISION array, dimension (min(M,N))
122*> The singular values of A, sorted so that S(i) >= S(i+1).
123*> \endverbatim
124*>
125*> \param[out] U
126*> \verbatim
127*> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
128*> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
129*> UCOL = min(M,N) if JOBZ = 'S'.
130*> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
131*> orthogonal matrix U;
132*> if JOBZ = 'S', U contains the first min(M,N) columns of U
133*> (the left singular vectors, stored columnwise);
134*> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
135*> \endverbatim
136*>
137*> \param[in] LDU
138*> \verbatim
139*> LDU is INTEGER
140*> The leading dimension of the array U. LDU >= 1; if
141*> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
142*> \endverbatim
143*>
144*> \param[out] VT
145*> \verbatim
146*> VT is DOUBLE PRECISION array, dimension (LDVT,N)
147*> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
148*> N-by-N orthogonal matrix V**T;
149*> if JOBZ = 'S', VT contains the first min(M,N) rows of
150*> V**T (the right singular vectors, stored rowwise);
151*> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
152*> \endverbatim
153*>
154*> \param[in] LDVT
155*> \verbatim
156*> LDVT is INTEGER
157*> The leading dimension of the array VT. LDVT >= 1;
158*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
159*> if JOBZ = 'S', LDVT >= min(M,N).
160*> \endverbatim
161*>
162*> \param[out] WORK
163*> \verbatim
164*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
166*> \endverbatim
167*>
168*> \param[in] LWORK
169*> \verbatim
170*> LWORK is INTEGER
171*> The dimension of the array WORK. LWORK >= 1.
172*> If LWORK = -1, a workspace query is assumed. The optimal
173*> size for the WORK array is calculated and stored in WORK(1),
174*> and no other work except argument checking is performed.
175*>
176*> Let mx = max(M,N) and mn = min(M,N).
177*> If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
178*> If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
179*> If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
180*> If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
181*> These are not tight minimums in all cases; see comments inside code.
182*> For good performance, LWORK should generally be larger;
183*> a query is recommended.
184*> \endverbatim
185*>
186*> \param[out] IWORK
187*> \verbatim
188*> IWORK is INTEGER array, dimension (8*min(M,N))
189*> \endverbatim
190*>
191*> \param[out] INFO
192*> \verbatim
193*> INFO is INTEGER
194*> < 0: if INFO = -i, the i-th argument had an illegal value.
195*> = -4: if A had a NAN entry.
196*> > 0: DBDSDC did not converge, updating process failed.
197*> = 0: successful exit.
198*> \endverbatim
199*
200* Authors:
201* ========
202*
203*> \author Univ. of Tennessee
204*> \author Univ. of California Berkeley
205*> \author Univ. of Colorado Denver
206*> \author NAG Ltd.
207*
208*> \ingroup doubleGEsing
209*
210*> \par Contributors:
211* ==================
212*>
213*> Ming Gu and Huan Ren, Computer Science Division, University of
214*> California at Berkeley, USA
215*>
216* =====================================================================
217 SUBROUTINE dgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
218 $ WORK, LWORK, IWORK, INFO )
219 implicit none
220*
221* -- LAPACK driver routine --
222* -- LAPACK is a software package provided by Univ. of Tennessee, --
223* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224*
225* .. Scalar Arguments ..
226 CHARACTER JOBZ
227 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
228* ..
229* .. Array Arguments ..
230 INTEGER IWORK( * )
231 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
232 $ vt( ldvt, * ), work( * )
233* ..
234*
235* =====================================================================
236*
237* .. Parameters ..
238 DOUBLE PRECISION ZERO, ONE
239 parameter( zero = 0.0d0, one = 1.0d0 )
240* ..
241* .. Local Scalars ..
242 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
243 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
244 $ ir, iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
245 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
246 $ mnthr, nwork, wrkbl
247 INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
248 $ lwork_dgebrd_nn, lwork_dgelqf_mn,
249 $ lwork_dgeqrf_mn,
250 $ lwork_dorgbr_p_mm, lwork_dorgbr_q_nn,
251 $ lwork_dorglq_mn, lwork_dorglq_nn,
252 $ lwork_dorgqr_mm, lwork_dorgqr_mn,
253 $ lwork_dormbr_prt_mm, lwork_dormbr_qln_mm,
254 $ lwork_dormbr_prt_mn, lwork_dormbr_qln_mn,
255 $ lwork_dormbr_prt_nn, lwork_dormbr_qln_nn
256 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
257* ..
258* .. Local Arrays ..
259 INTEGER IDUM( 1 )
260 DOUBLE PRECISION DUM( 1 )
261* ..
262* .. External Subroutines ..
263 EXTERNAL dbdsdc, dgebrd, dgelqf, dgemm, dgeqrf, dlacpy,
264 $ dlascl, dlaset, dorgbr, dorglq, dorgqr, dormbr,
265 $ xerbla
266* ..
267* .. External Functions ..
268 LOGICAL LSAME, DISNAN
269 DOUBLE PRECISION DLAMCH, DLANGE, DROUNDUP_LWORK
270 EXTERNAL dlamch, dlange, lsame, disnan,
271 $ droundup_lwork
272* ..
273* .. Intrinsic Functions ..
274 INTRINSIC int, max, min, sqrt
275* ..
276* .. Executable Statements ..
277*
278* Test the input arguments
279*
280 info = 0
281 minmn = min( m, n )
282 wntqa = lsame( jobz, 'A' )
283 wntqs = lsame( jobz, 'S' )
284 wntqas = wntqa .OR. wntqs
285 wntqo = lsame( jobz, 'O' )
286 wntqn = lsame( jobz, 'N' )
287 lquery = ( lwork.EQ.-1 )
288*
289 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) ) THEN
290 info = -1
291 ELSE IF( m.LT.0 ) THEN
292 info = -2
293 ELSE IF( n.LT.0 ) THEN
294 info = -3
295 ELSE IF( lda.LT.max( 1, m ) ) THEN
296 info = -5
297 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
298 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) ) THEN
299 info = -8
300 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
301 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
302 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) ) THEN
303 info = -10
304 END IF
305*
306* Compute workspace
307* Note: Comments in the code beginning "Workspace:" describe the
308* minimal amount of workspace allocated at that point in the code,
309* as well as the preferred amount for good performance.
310* NB refers to the optimal block size for the immediately
311* following subroutine, as returned by ILAENV.
312*
313 IF( info.EQ.0 ) THEN
314 minwrk = 1
315 maxwrk = 1
316 bdspac = 0
317 mnthr = int( minmn*11.0d0 / 6.0d0 )
318 IF( m.GE.n .AND. minmn.GT.0 ) THEN
319*
320* Compute space needed for DBDSDC
321*
322 IF( wntqn ) THEN
323* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
324* keep 7*N for backwards compatibility.
325 bdspac = 7*n
326 ELSE
327 bdspac = 3*n*n + 4*n
328 END IF
329*
330* Compute space preferred for each routine
331 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
332 $ dum(1), dum(1), -1, ierr )
333 lwork_dgebrd_mn = int( dum(1) )
334*
335 CALL dgebrd( n, n, dum(1), n, dum(1), dum(1), dum(1),
336 $ dum(1), dum(1), -1, ierr )
337 lwork_dgebrd_nn = int( dum(1) )
338*
339 CALL dgeqrf( m, n, dum(1), m, dum(1), dum(1), -1, ierr )
340 lwork_dgeqrf_mn = int( dum(1) )
341*
342 CALL dorgbr( 'Q', n, n, n, dum(1), n, dum(1), dum(1), -1,
343 $ ierr )
344 lwork_dorgbr_q_nn = int( dum(1) )
345*
346 CALL dorgqr( m, m, n, dum(1), m, dum(1), dum(1), -1, ierr )
347 lwork_dorgqr_mm = int( dum(1) )
348*
349 CALL dorgqr( m, n, n, dum(1), m, dum(1), dum(1), -1, ierr )
350 lwork_dorgqr_mn = int( dum(1) )
351*
352 CALL dormbr( 'P', 'R', 'T', n, n, n, dum(1), n,
353 $ dum(1), dum(1), n, dum(1), -1, ierr )
354 lwork_dormbr_prt_nn = int( dum(1) )
355*
356 CALL dormbr( 'Q', 'L', 'N', n, n, n, dum(1), n,
357 $ dum(1), dum(1), n, dum(1), -1, ierr )
358 lwork_dormbr_qln_nn = int( dum(1) )
359*
360 CALL dormbr( 'Q', 'L', 'N', m, n, n, dum(1), m,
361 $ dum(1), dum(1), m, dum(1), -1, ierr )
362 lwork_dormbr_qln_mn = int( dum(1) )
363*
364 CALL dormbr( 'Q', 'L', 'N', m, m, n, dum(1), m,
365 $ dum(1), dum(1), m, dum(1), -1, ierr )
366 lwork_dormbr_qln_mm = int( dum(1) )
367*
368 IF( m.GE.mnthr ) THEN
369 IF( wntqn ) THEN
370*
371* Path 1 (M >> N, JOBZ='N')
372*
373 wrkbl = n + lwork_dgeqrf_mn
374 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
375 maxwrk = max( wrkbl, bdspac + n )
376 minwrk = bdspac + n
377 ELSE IF( wntqo ) THEN
378*
379* Path 2 (M >> N, JOBZ='O')
380*
381 wrkbl = n + lwork_dgeqrf_mn
382 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
383 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
384 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
385 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
386 wrkbl = max( wrkbl, 3*n + bdspac )
387 maxwrk = wrkbl + 2*n*n
388 minwrk = bdspac + 2*n*n + 3*n
389 ELSE IF( wntqs ) THEN
390*
391* Path 3 (M >> N, JOBZ='S')
392*
393 wrkbl = n + lwork_dgeqrf_mn
394 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
395 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
396 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
397 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
398 wrkbl = max( wrkbl, 3*n + bdspac )
399 maxwrk = wrkbl + n*n
400 minwrk = bdspac + n*n + 3*n
401 ELSE IF( wntqa ) THEN
402*
403* Path 4 (M >> N, JOBZ='A')
404*
405 wrkbl = n + lwork_dgeqrf_mn
406 wrkbl = max( wrkbl, n + lwork_dorgqr_mm )
407 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
408 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
409 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
410 wrkbl = max( wrkbl, 3*n + bdspac )
411 maxwrk = wrkbl + n*n
412 minwrk = n*n + max( 3*n + bdspac, n + m )
413 END IF
414 ELSE
415*
416* Path 5 (M >= N, but not much larger)
417*
418 wrkbl = 3*n + lwork_dgebrd_mn
419 IF( wntqn ) THEN
420* Path 5n (M >= N, jobz='N')
421 maxwrk = max( wrkbl, 3*n + bdspac )
422 minwrk = 3*n + max( m, bdspac )
423 ELSE IF( wntqo ) THEN
424* Path 5o (M >= N, jobz='O')
425 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
426 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
427 wrkbl = max( wrkbl, 3*n + bdspac )
428 maxwrk = wrkbl + m*n
429 minwrk = 3*n + max( m, n*n + bdspac )
430 ELSE IF( wntqs ) THEN
431* Path 5s (M >= N, jobz='S')
432 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
433 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
434 maxwrk = max( wrkbl, 3*n + bdspac )
435 minwrk = 3*n + max( m, bdspac )
436 ELSE IF( wntqa ) THEN
437* Path 5a (M >= N, jobz='A')
438 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mm )
439 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
440 maxwrk = max( wrkbl, 3*n + bdspac )
441 minwrk = 3*n + max( m, bdspac )
442 END IF
443 END IF
444 ELSE IF( minmn.GT.0 ) THEN
445*
446* Compute space needed for DBDSDC
447*
448 IF( wntqn ) THEN
449* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
450* keep 7*N for backwards compatibility.
451 bdspac = 7*m
452 ELSE
453 bdspac = 3*m*m + 4*m
454 END IF
455*
456* Compute space preferred for each routine
457 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
458 $ dum(1), dum(1), -1, ierr )
459 lwork_dgebrd_mn = int( dum(1) )
460*
461 CALL dgebrd( m, m, a, m, s, dum(1), dum(1),
462 $ dum(1), dum(1), -1, ierr )
463 lwork_dgebrd_mm = int( dum(1) )
464*
465 CALL dgelqf( m, n, a, m, dum(1), dum(1), -1, ierr )
466 lwork_dgelqf_mn = int( dum(1) )
467*
468 CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
469 lwork_dorglq_nn = int( dum(1) )
470*
471 CALL dorglq( m, n, m, a, m, dum(1), dum(1), -1, ierr )
472 lwork_dorglq_mn = int( dum(1) )
473*
474 CALL dorgbr( 'P', m, m, m, a, n, dum(1), dum(1), -1, ierr )
475 lwork_dorgbr_p_mm = int( dum(1) )
476*
477 CALL dormbr( 'P', 'R', 'T', m, m, m, dum(1), m,
478 $ dum(1), dum(1), m, dum(1), -1, ierr )
479 lwork_dormbr_prt_mm = int( dum(1) )
480*
481 CALL dormbr( 'P', 'R', 'T', m, n, m, dum(1), m,
482 $ dum(1), dum(1), m, dum(1), -1, ierr )
483 lwork_dormbr_prt_mn = int( dum(1) )
484*
485 CALL dormbr( 'P', 'R', 'T', n, n, m, dum(1), n,
486 $ dum(1), dum(1), n, dum(1), -1, ierr )
487 lwork_dormbr_prt_nn = int( dum(1) )
488*
489 CALL dormbr( 'Q', 'L', 'N', m, m, m, dum(1), m,
490 $ dum(1), dum(1), m, dum(1), -1, ierr )
491 lwork_dormbr_qln_mm = int( dum(1) )
492*
493 IF( n.GE.mnthr ) THEN
494 IF( wntqn ) THEN
495*
496* Path 1t (N >> M, JOBZ='N')
497*
498 wrkbl = m + lwork_dgelqf_mn
499 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
500 maxwrk = max( wrkbl, bdspac + m )
501 minwrk = bdspac + m
502 ELSE IF( wntqo ) THEN
503*
504* Path 2t (N >> M, JOBZ='O')
505*
506 wrkbl = m + lwork_dgelqf_mn
507 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
508 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
509 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
510 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
511 wrkbl = max( wrkbl, 3*m + bdspac )
512 maxwrk = wrkbl + 2*m*m
513 minwrk = bdspac + 2*m*m + 3*m
514 ELSE IF( wntqs ) THEN
515*
516* Path 3t (N >> M, JOBZ='S')
517*
518 wrkbl = m + lwork_dgelqf_mn
519 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
520 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
521 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
522 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
523 wrkbl = max( wrkbl, 3*m + bdspac )
524 maxwrk = wrkbl + m*m
525 minwrk = bdspac + m*m + 3*m
526 ELSE IF( wntqa ) THEN
527*
528* Path 4t (N >> M, JOBZ='A')
529*
530 wrkbl = m + lwork_dgelqf_mn
531 wrkbl = max( wrkbl, m + lwork_dorglq_nn )
532 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
533 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
534 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
535 wrkbl = max( wrkbl, 3*m + bdspac )
536 maxwrk = wrkbl + m*m
537 minwrk = m*m + max( 3*m + bdspac, m + n )
538 END IF
539 ELSE
540*
541* Path 5t (N > M, but not much larger)
542*
543 wrkbl = 3*m + lwork_dgebrd_mn
544 IF( wntqn ) THEN
545* Path 5tn (N > M, jobz='N')
546 maxwrk = max( wrkbl, 3*m + bdspac )
547 minwrk = 3*m + max( n, bdspac )
548 ELSE IF( wntqo ) THEN
549* Path 5to (N > M, jobz='O')
550 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
551 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
552 wrkbl = max( wrkbl, 3*m + bdspac )
553 maxwrk = wrkbl + m*n
554 minwrk = 3*m + max( n, m*m + bdspac )
555 ELSE IF( wntqs ) THEN
556* Path 5ts (N > M, jobz='S')
557 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
558 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
559 maxwrk = max( wrkbl, 3*m + bdspac )
560 minwrk = 3*m + max( n, bdspac )
561 ELSE IF( wntqa ) THEN
562* Path 5ta (N > M, jobz='A')
563 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
564 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_nn )
565 maxwrk = max( wrkbl, 3*m + bdspac )
566 minwrk = 3*m + max( n, bdspac )
567 END IF
568 END IF
569 END IF
570
571 maxwrk = max( maxwrk, minwrk )
572 work( 1 ) = droundup_lwork( maxwrk )
573*
574 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
575 info = -12
576 END IF
577 END IF
578*
579 IF( info.NE.0 ) THEN
580 CALL xerbla( 'DGESDD', -info )
581 RETURN
582 ELSE IF( lquery ) THEN
583 RETURN
584 END IF
585*
586* Quick return if possible
587*
588 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
589 RETURN
590 END IF
591*
592* Get machine constants
593*
594 eps = dlamch( 'P' )
595 smlnum = sqrt( dlamch( 'S' ) ) / eps
596 bignum = one / smlnum
597*
598* Scale A if max element outside range [SMLNUM,BIGNUM]
599*
600 anrm = dlange( 'M', m, n, a, lda, dum )
601 IF( disnan( anrm ) ) THEN
602 info = -4
603 RETURN
604 END IF
605 iscl = 0
606 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
607 iscl = 1
608 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
609 ELSE IF( anrm.GT.bignum ) THEN
610 iscl = 1
611 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
612 END IF
613*
614 IF( m.GE.n ) THEN
615*
616* A has at least as many rows as columns. If A has sufficiently
617* more rows than columns, first reduce using the QR
618* decomposition (if sufficient workspace available)
619*
620 IF( m.GE.mnthr ) THEN
621*
622 IF( wntqn ) THEN
623*
624* Path 1 (M >> N, JOBZ='N')
625* No singular vectors to be computed
626*
627 itau = 1
628 nwork = itau + n
629*
630* Compute A=Q*R
631* Workspace: need N [tau] + N [work]
632* Workspace: prefer N [tau] + N*NB [work]
633*
634 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
635 $ lwork - nwork + 1, ierr )
636*
637* Zero out below R
638*
639 CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
640 ie = 1
641 itauq = ie + n
642 itaup = itauq + n
643 nwork = itaup + n
644*
645* Bidiagonalize R in A
646* Workspace: need 3*N [e, tauq, taup] + N [work]
647* Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
648*
649 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
650 $ work( itaup ), work( nwork ), lwork-nwork+1,
651 $ ierr )
652 nwork = ie + n
653*
654* Perform bidiagonal SVD, computing singular values only
655* Workspace: need N [e] + BDSPAC
656*
657 CALL dbdsdc( 'U', 'N', n, s, work( ie ), dum, 1, dum, 1,
658 $ dum, idum, work( nwork ), iwork, info )
659*
660 ELSE IF( wntqo ) THEN
661*
662* Path 2 (M >> N, JOBZ = 'O')
663* N left singular vectors to be overwritten on A and
664* N right singular vectors to be computed in VT
665*
666 ir = 1
667*
668* WORK(IR) is LDWRKR by N
669*
670 IF( lwork .GE. lda*n + n*n + 3*n + bdspac ) THEN
671 ldwrkr = lda
672 ELSE
673 ldwrkr = ( lwork - n*n - 3*n - bdspac ) / n
674 END IF
675 itau = ir + ldwrkr*n
676 nwork = itau + n
677*
678* Compute A=Q*R
679* Workspace: need N*N [R] + N [tau] + N [work]
680* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
681*
682 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
683 $ lwork - nwork + 1, ierr )
684*
685* Copy R to WORK(IR), zeroing out below it
686*
687 CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
688 CALL dlaset( 'L', n - 1, n - 1, zero, zero, work(ir+1),
689 $ ldwrkr )
690*
691* Generate Q in A
692* Workspace: need N*N [R] + N [tau] + N [work]
693* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
694*
695 CALL dorgqr( m, n, n, a, lda, work( itau ),
696 $ work( nwork ), lwork - nwork + 1, ierr )
697 ie = itau
698 itauq = ie + n
699 itaup = itauq + n
700 nwork = itaup + n
701*
702* Bidiagonalize R in WORK(IR)
703* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
704* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
705*
706 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
707 $ work( itauq ), work( itaup ), work( nwork ),
708 $ lwork - nwork + 1, ierr )
709*
710* WORK(IU) is N by N
711*
712 iu = nwork
713 nwork = iu + n*n
714*
715* Perform bidiagonal SVD, computing left singular vectors
716* of bidiagonal matrix in WORK(IU) and computing right
717* singular vectors of bidiagonal matrix in VT
718* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
719*
720 CALL dbdsdc( 'U', 'I', n, s, work( ie ), work( iu ), n,
721 $ vt, ldvt, dum, idum, work( nwork ), iwork,
722 $ info )
723*
724* Overwrite WORK(IU) by left singular vectors of R
725* and VT by right singular vectors of R
726* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N [work]
727* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
728*
729 CALL dormbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
730 $ work( itauq ), work( iu ), n, work( nwork ),
731 $ lwork - nwork + 1, ierr )
732 CALL dormbr( 'P', 'R', 'T', n, n, n, work( ir ), ldwrkr,
733 $ work( itaup ), vt, ldvt, work( nwork ),
734 $ lwork - nwork + 1, ierr )
735*
736* Multiply Q in A by left singular vectors of R in
737* WORK(IU), storing result in WORK(IR) and copying to A
738* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U]
739* Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U]
740*
741 DO 10 i = 1, m, ldwrkr
742 chunk = min( m - i + 1, ldwrkr )
743 CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
744 $ lda, work( iu ), n, zero, work( ir ),
745 $ ldwrkr )
746 CALL dlacpy( 'F', chunk, n, work( ir ), ldwrkr,
747 $ a( i, 1 ), lda )
748 10 CONTINUE
749*
750 ELSE IF( wntqs ) THEN
751*
752* Path 3 (M >> N, JOBZ='S')
753* N left singular vectors to be computed in U and
754* N right singular vectors to be computed in VT
755*
756 ir = 1
757*
758* WORK(IR) is N by N
759*
760 ldwrkr = n
761 itau = ir + ldwrkr*n
762 nwork = itau + n
763*
764* Compute A=Q*R
765* Workspace: need N*N [R] + N [tau] + N [work]
766* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
767*
768 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
769 $ lwork - nwork + 1, ierr )
770*
771* Copy R to WORK(IR), zeroing out below it
772*
773 CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
774 CALL dlaset( 'L', n - 1, n - 1, zero, zero, work(ir+1),
775 $ ldwrkr )
776*
777* Generate Q in A
778* Workspace: need N*N [R] + N [tau] + N [work]
779* Workspace: prefer N*N [R] + N [tau] + N*NB [work]
780*
781 CALL dorgqr( m, n, n, a, lda, work( itau ),
782 $ work( nwork ), lwork - nwork + 1, ierr )
783 ie = itau
784 itauq = ie + n
785 itaup = itauq + n
786 nwork = itaup + n
787*
788* Bidiagonalize R in WORK(IR)
789* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
790* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
791*
792 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
793 $ work( itauq ), work( itaup ), work( nwork ),
794 $ lwork - nwork + 1, ierr )
795*
796* Perform bidiagonal SVD, computing left singular vectors
797* of bidiagoal matrix in U and computing right singular
798* vectors of bidiagonal matrix in VT
799* Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC
800*
801 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
802 $ ldvt, dum, idum, work( nwork ), iwork,
803 $ info )
804*
805* Overwrite U by left singular vectors of R and VT
806* by right singular vectors of R
807* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work]
808* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work]
809*
810 CALL dormbr( 'Q', 'l', 'n', N, N, N, WORK( IR ), LDWRKR,
811 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
812 $ LWORK - NWORK + 1, IERR )
813*
814 CALL DORMBR( 'p', 'r', 't', N, N, N, WORK( IR ), LDWRKR,
815 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
816 $ LWORK - NWORK + 1, IERR )
817*
818* Multiply Q in A by left singular vectors of R in
819* WORK(IR), storing result in U
820* Workspace: need N*N [R]
821*
822 CALL DLACPY( 'f', N, N, U, LDU, WORK( IR ), LDWRKR )
823 CALL DGEMM( 'n', 'n', M, N, N, ONE, A, LDA, WORK( IR ),
824 $ LDWRKR, ZERO, U, LDU )
825*
826 ELSE IF( WNTQA ) THEN
827*
828* Path 4 (M >> N, JOBZ='A')
829* M left singular vectors to be computed in U and
830* N right singular vectors to be computed in VT
831*
832 IU = 1
833*
834* WORK(IU) is N by N
835*
836 LDWRKU = N
837 ITAU = IU + LDWRKU*N
838 NWORK = ITAU + N
839*
840* Compute A=Q*R, copying result to U
841* Workspace: need N*N [U] + N [tau] + N [work]
842* Workspace: prefer N*N [U] + N [tau] + N*NB [work]
843*
844 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
845 $ LWORK - NWORK + 1, IERR )
846 CALL DLACPY( 'l', M, N, A, LDA, U, LDU )
847*
848* Generate Q in U
849* Workspace: need N*N [U] + N [tau] + M [work]
850* Workspace: prefer N*N [U] + N [tau] + M*NB [work]
851 CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
852 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
853*
854* Produce R in A, zeroing out other entries
855*
856 CALL DLASET( 'l', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
857 IE = ITAU
858 ITAUQ = IE + N
859 ITAUP = ITAUQ + N
860 NWORK = ITAUP + N
861*
862* Bidiagonalize R in A
863* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
864* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work]
865*
866 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
867 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
868 $ IERR )
869*
870* Perform bidiagonal SVD, computing left singular vectors
871* of bidiagonal matrix in WORK(IU) and computing right
872* singular vectors of bidiagonal matrix in VT
873* Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC
874*
875 CALL DBDSDC( 'u', 'i', N, S, WORK( IE ), WORK( IU ), N,
876 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
877 $ INFO )
878*
879* Overwrite WORK(IU) by left singular vectors of R and VT
880* by right singular vectors of R
881* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work]
882* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work]
883*
884 CALL DORMBR( 'q', 'l', 'n', N, N, N, A, LDA,
885 $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
886 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
887 CALL DORMBR( 'p', 'r', 't', N, N, N, A, LDA,
888 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
889 $ LWORK - NWORK + 1, IERR )
890*
891* Multiply Q in U by left singular vectors of R in
892* WORK(IU), storing result in A
893* Workspace: need N*N [U]
894*
895 CALL DGEMM( 'n', 'n', M, N, N, ONE, U, LDU, WORK( IU ),
896 $ LDWRKU, ZERO, A, LDA )
897*
898* Copy left singular vectors of A from A to U
899*
900 CALL DLACPY( 'f', M, N, A, LDA, U, LDU )
901*
902 END IF
903*
904 ELSE
905*
906* M .LT. MNTHR
907*
908* Path 5 (M >= N, but not much larger)
909* Reduce to bidiagonal form without QR decomposition
910*
911 IE = 1
912 ITAUQ = IE + N
913 ITAUP = ITAUQ + N
914 NWORK = ITAUP + N
915*
916* Bidiagonalize A
917* Workspace: need 3*N [e, tauq, taup] + M [work]
918* Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
919*
920 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
921 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
922 $ IERR )
923 IF( WNTQN ) THEN
924*
925* Path 5n (M >= N, JOBZ='N')
926* Perform bidiagonal SVD, only computing singular values
927* Workspace: need 3*N [e, tauq, taup] + BDSPAC
928*
929 CALL DBDSDC( 'u', 'n', N, S, WORK( IE ), DUM, 1, DUM, 1,
930 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
931 ELSE IF( WNTQO ) THEN
932* Path 5o (M >= N, JOBZ='O')
933 IU = NWORK
934.GE. IF( LWORK M*N + 3*N + BDSPAC ) THEN
935*
936* WORK( IU ) is M by N
937*
938 LDWRKU = M
939 NWORK = IU + LDWRKU*N
940 CALL DLASET( 'f', M, N, ZERO, ZERO, WORK( IU ),
941 $ LDWRKU )
942* IR is unused; silence compile warnings
943 IR = -1
944 ELSE
945*
946* WORK( IU ) is N by N
947*
948 LDWRKU = N
949 NWORK = IU + LDWRKU*N
950*
951* WORK(IR) is LDWRKR by N
952*
953 IR = NWORK
954 LDWRKR = ( LWORK - N*N - 3*N ) / N
955 END IF
956 NWORK = IU + LDWRKU*N
957*
958* Perform bidiagonal SVD, computing left singular vectors
959* of bidiagonal matrix in WORK(IU) and computing right
960* singular vectors of bidiagonal matrix in VT
961* Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC
962*
963 CALL DBDSDC( 'u', 'i', N, S, WORK( IE ), WORK( IU ),
964 $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
965 $ IWORK, INFO )
966*
967* Overwrite VT by right singular vectors of A
968* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
969* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
970*
971 CALL DORMBR( 'p', 'r', 't', N, N, N, A, LDA,
972 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
973 $ LWORK - NWORK + 1, IERR )
974*
975.GE. IF( LWORK M*N + 3*N + BDSPAC ) THEN
976*
977* Path 5o-fast
978* Overwrite WORK(IU) by left singular vectors of A
979* Workspace: need 3*N [e, tauq, taup] + M*N [U] + N [work]
980* Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work]
981*
982 CALL DORMBR( 'q', 'l', 'N', m, n, n, a, lda,
983 $ work( itauq ), work( iu ), ldwrku,
984 $ work( nwork ), lwork - nwork + 1, ierr )
985*
986* Copy left singular vectors of A from WORK(IU) to A
987*
988 CALL dlacpy( 'F', m, n, work( iu ), ldwrku, a, lda )
989 ELSE
990*
991* Path 5o-slow
992* Generate Q in A
993* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work]
994* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
995*
996 CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
997 $ work( nwork ), lwork - nwork + 1, ierr )
998*
999* Multiply Q in A by left singular vectors of
1000* bidiagonal matrix in WORK(IU), storing result in
1001* WORK(IR) and copying to A
1002* Workspace: need 3*N [e, tauq, taup] + N*N [U] + NB*N [R]
1003* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N [R]
1004*
1005 DO 20 i = 1, m, ldwrkr
1006 chunk = min( m - i + 1, ldwrkr )
1007 CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
1008 $ lda, work( iu ), ldwrku, zero,
1009 $ work( ir ), ldwrkr )
1010 CALL dlacpy( 'F', chunk, n, work( ir ), ldwrkr,
1011 $ a( i, 1 ), lda )
1012 20 CONTINUE
1013 END IF
1014*
1015 ELSE IF( wntqs ) THEN
1016*
1017* Path 5s (M >= N, JOBZ='S')
1018* Perform bidiagonal SVD, computing left singular vectors
1019* of bidiagonal matrix in U and computing right singular
1020* vectors of bidiagonal matrix in VT
1021* Workspace: need 3*N [e, tauq, taup] + BDSPAC
1022*
1023 CALL dlaset( 'F', m, n, zero, zero, u, ldu )
1024 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
1025 $ ldvt, dum, idum, work( nwork ), iwork,
1026 $ info )
1027*
1028* Overwrite U by left singular vectors of A and VT
1029* by right singular vectors of A
1030* Workspace: need 3*N [e, tauq, taup] + N [work]
1031* Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1032*
1033 CALL dormbr( 'Q', 'L', 'N', m, n, n, a, lda,
1034 $ work( itauq ), u, ldu, work( nwork ),
1035 $ lwork - nwork + 1, ierr )
1036 CALL dormbr( 'P', 'R', 'T', n, n, n, a, lda,
1037 $ work( itaup ), vt, ldvt, work( nwork ),
1038 $ lwork - nwork + 1, ierr )
1039 ELSE IF( wntqa ) THEN
1040*
1041* Path 5a (M >= N, JOBZ='A')
1042* Perform bidiagonal SVD, computing left singular vectors
1043* of bidiagonal matrix in U and computing right singular
1044* vectors of bidiagonal matrix in VT
1045* Workspace: need 3*N [e, tauq, taup] + BDSPAC
1046*
1047 CALL dlaset( 'F', m, m, zero, zero, u, ldu )
1048 CALL dbdsdc( 'U', 'I', n, s, work( ie ), u, ldu, vt,
1049 $ ldvt, dum, idum, work( nwork ), iwork,
1050 $ info )
1051*
1052* Set the right corner of U to identity matrix
1053*
1054 IF( m.GT.n ) THEN
1055 CALL dlaset( 'F', m - n, m - n, zero, one, u(n+1,n+1),
1056 $ ldu )
1057 END IF
1058*
1059* Overwrite U by left singular vectors of A and VT
1060* by right singular vectors of A
1061* Workspace: need 3*N [e, tauq, taup] + M [work]
1062* Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1063*
1064 CALL dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1065 $ work( itauq ), u, ldu, work( nwork ),
1066 $ lwork - nwork + 1, ierr )
1067 CALL dormbr( 'P', 'R', 'T', n, n, m, a, lda,
1068 $ work( itaup ), vt, ldvt, work( nwork ),
1069 $ lwork - nwork + 1, ierr )
1070 END IF
1071*
1072 END IF
1073*
1074 ELSE
1075*
1076* A has more columns than rows. If A has sufficiently more
1077* columns than rows, first reduce using the LQ decomposition (if
1078* sufficient workspace available)
1079*
1080 IF( n.GE.mnthr ) THEN
1081*
1082 IF( wntqn ) THEN
1083*
1084* Path 1t (N >> M, JOBZ='N')
1085* No singular vectors to be computed
1086*
1087 itau = 1
1088 nwork = itau + m
1089*
1090* Compute A=L*Q
1091* Workspace: need M [tau] + M [work]
1092* Workspace: prefer M [tau] + M*NB [work]
1093*
1094 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1095 $ lwork - nwork + 1, ierr )
1096*
1097* Zero out above L
1098*
1099 CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1100 ie = 1
1101 itauq = ie + m
1102 itaup = itauq + m
1103 nwork = itaup + m
1104*
1105* Bidiagonalize L in A
1106* Workspace: need 3*M [e, tauq, taup] + M [work]
1107* Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1108*
1109 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1110 $ work( itaup ), work( nwork ), lwork-nwork+1,
1111 $ ierr )
1112 nwork = ie + m
1113*
1114* Perform bidiagonal SVD, computing singular values only
1115* Workspace: need M [e] + BDSPAC
1116*
1117 CALL dbdsdc( 'U', 'N', m, s, work( ie ), dum, 1, dum, 1,
1118 $ dum, idum, work( nwork ), iwork, info )
1119*
1120 ELSE IF( wntqo ) THEN
1121*
1122* Path 2t (N >> M, JOBZ='O')
1123* M right singular vectors to be overwritten on A and
1124* M left singular vectors to be computed in U
1125*
1126 ivt = 1
1127*
1128* WORK(IVT) is M by M
1129* WORK(IL) is M by M; it is later resized to M by chunk for gemm
1130*
1131 il = ivt + m*m
1132 IF( lwork .GE. m*n + m*m + 3*m + bdspac ) THEN
1133 ldwrkl = m
1134 chunk = n
1135 ELSE
1136 ldwrkl = m
1137 chunk = ( lwork - m*m ) / m
1138 END IF
1139 itau = il + ldwrkl*m
1140 nwork = itau + m
1141*
1142* Compute A=L*Q
1143* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1144* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1145*
1146 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1147 $ lwork - nwork + 1, ierr )
1148*
1149* Copy L to WORK(IL), zeroing about above it
1150*
1151 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1152 CALL dlaset( 'U', m - 1, m - 1, zero, zero,
1153 $ work( il + ldwrkl ), ldwrkl )
1154*
1155* Generate Q in A
1156* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1157* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1158*
1159 CALL dorglq( m, n, m, a, lda, work( itau ),
1160 $ work( nwork ), lwork - nwork + 1, ierr )
1161 ie = itau
1162 itauq = ie + m
1163 itaup = itauq + m
1164 nwork = itaup + m
1165*
1166* Bidiagonalize L in WORK(IL)
1167* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1168* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1169*
1170 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1171 $ work( itauq ), work( itaup ), work( nwork ),
1172 $ lwork - nwork + 1, ierr )
1173*
1174* Perform bidiagonal SVD, computing left singular vectors
1175* of bidiagonal matrix in U, and computing right singular
1176* vectors of bidiagonal matrix in WORK(IVT)
1177* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1178*
1179 CALL dbdsdc( 'U', 'I', m, s, work( ie ), u, ldu,
1180 $ work( ivt ), m, dum, idum, work( nwork ),
1181 $ iwork, info )
1182*
1183* Overwrite U by left singular vectors of L and WORK(IVT)
1184* by right singular vectors of L
1185* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work]
1186* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1187*
1188 CALL dormbr( 'Q', 'l', 'n', M, M, M, WORK( IL ), LDWRKL,
1189 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1190 $ LWORK - NWORK + 1, IERR )
1191 CALL DORMBR( 'p', 'r', 't', M, M, M, WORK( IL ), LDWRKL,
1192 $ WORK( ITAUP ), WORK( IVT ), M,
1193 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1194*
1195* Multiply right singular vectors of L in WORK(IVT) by Q
1196* in A, storing result in WORK(IL) and copying to A
1197* Workspace: need M*M [VT] + M*M [L]
1198* Workspace: prefer M*M [VT] + M*N [L]
1199* At this point, L is resized as M by chunk.
1200*
1201 DO 30 I = 1, N, CHUNK
1202 BLK = MIN( N - I + 1, CHUNK )
1203 CALL DGEMM( 'n', 'n', M, BLK, M, ONE, WORK( IVT ), M,
1204 $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1205 CALL DLACPY( 'f', M, BLK, WORK( IL ), LDWRKL,
1206 $ A( 1, I ), LDA )
1207 30 CONTINUE
1208*
1209 ELSE IF( WNTQS ) THEN
1210*
1211* Path 3t (N >> M, JOBZ='S')
1212* M right singular vectors to be computed in VT and
1213* M left singular vectors to be computed in U
1214*
1215 IL = 1
1216*
1217* WORK(IL) is M by M
1218*
1219 LDWRKL = M
1220 ITAU = IL + LDWRKL*M
1221 NWORK = ITAU + M
1222*
1223* Compute A=L*Q
1224* Workspace: need M*M [L] + M [tau] + M [work]
1225* Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1226*
1227 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1228 $ LWORK - NWORK + 1, IERR )
1229*
1230* Copy L to WORK(IL), zeroing out above it
1231*
1232 CALL DLACPY( 'l', M, M, A, LDA, WORK( IL ), LDWRKL )
1233 CALL DLASET( 'u', M - 1, M - 1, ZERO, ZERO,
1234 $ WORK( IL + LDWRKL ), LDWRKL )
1235*
1236* Generate Q in A
1237* Workspace: need M*M [L] + M [tau] + M [work]
1238* Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1239*
1240 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1241 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1242 IE = ITAU
1243 ITAUQ = IE + M
1244 ITAUP = ITAUQ + M
1245 NWORK = ITAUP + M
1246*
1247* Bidiagonalize L in WORK(IU).
1248* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1249* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1250*
1251 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1252 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1253 $ LWORK - NWORK + 1, IERR )
1254*
1255* Perform bidiagonal SVD, computing left singular vectors
1256* of bidiagonal matrix in U and computing right singular
1257* vectors of bidiagonal matrix in VT
1258* Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1259*
1260 CALL DBDSDC( 'u', 'i', M, S, WORK( IE ), U, LDU, VT,
1261 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1262 $ INFO )
1263*
1264* Overwrite U by left singular vectors of L and VT
1265* by right singular vectors of L
1266* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work]
1267* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1268*
1269 CALL DORMBR( 'q', 'l', 'n', M, M, M, WORK( IL ), LDWRKL,
1270 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1271 $ LWORK - NWORK + 1, IERR )
1272 CALL DORMBR( 'p', 'r', 't', M, M, M, WORK( IL ), LDWRKL,
1273 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1274 $ LWORK - NWORK + 1, IERR )
1275*
1276* Multiply right singular vectors of L in WORK(IL) by
1277* Q in A, storing result in VT
1278* Workspace: need M*M [L]
1279*
1280 CALL DLACPY( 'f', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1281 CALL DGEMM( 'n', 'n', M, N, M, ONE, WORK( IL ), LDWRKL,
1282 $ A, LDA, ZERO, VT, LDVT )
1283*
1284 ELSE IF( WNTQA ) THEN
1285*
1286* Path 4t (N >> M, JOBZ='A')
1287* N right singular vectors to be computed in VT and
1288* M left singular vectors to be computed in U
1289*
1290 IVT = 1
1291*
1292* WORK(IVT) is M by M
1293*
1294 LDWKVT = M
1295 ITAU = IVT + LDWKVT*M
1296 NWORK = ITAU + M
1297*
1298* Compute A=L*Q, copying result to VT
1299* Workspace: need M*M [VT] + M [tau] + M [work]
1300* Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1301*
1302 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1303 $ LWORK - NWORK + 1, IERR )
1304 CALL DLACPY( 'u', M, N, A, LDA, VT, LDVT )
1305*
1306* Generate Q in VT
1307* Workspace: need M*M [VT] + M [tau] + N [work]
1308* Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1309*
1310 CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1311 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1312*
1313* Produce L in A, zeroing out other entries
1314*
1315 CALL DLASET( 'u', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1316 IE = ITAU
1317 ITAUQ = IE + M
1318 ITAUP = ITAUQ + M
1319 NWORK = ITAUP + M
1320*
1321* Bidiagonalize L in A
1322* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + M [work]
1323* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work]
1324*
1325 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1326 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1327 $ IERR )
1328*
1329* Perform bidiagonal SVD, computing left singular vectors
1330* of bidiagonal matrix in U and computing right singular
1331* vectors of bidiagonal matrix in WORK(IVT)
1332* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1333*
1334 CALL DBDSDC( 'u', 'i', M, S, WORK( IE ), U, LDU,
1335 $ WORK( IVT ), LDWKVT, DUM, IDUM,
1336 $ WORK( NWORK ), IWORK, INFO )
1337*
1338* Overwrite U by left singular vectors of L and WORK(IVT)
1339* by right singular vectors of L
1340* Workspace: need M*M [VT] + 3*M [e, tauq, taup]+ M [work]
1341* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work]
1342*
1343 CALL DORMBR( 'q', 'l', 'n', M, M, M, A, LDA,
1344 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1345 $ LWORK - NWORK + 1, IERR )
1346 CALL DORMBR( 'p', 'r', 't', M, M, M, A, LDA,
1347 $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1348 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1349*
1350* Multiply right singular vectors of L in WORK(IVT) by
1351* Q in VT, storing result in A
1352* Workspace: need M*M [VT]
1353*
1354 CALL DGEMM( 'n', 'n', M, N, M, ONE, WORK( IVT ), LDWKVT,
1355 $ VT, LDVT, ZERO, A, LDA )
1356*
1357* Copy right singular vectors of A from A to VT
1358*
1359 CALL DLACPY( 'f', M, N, A, LDA, VT, LDVT )
1360*
1361 END IF
1362*
1363 ELSE
1364*
1365* N .LT. MNTHR
1366*
1367* Path 5t (N > M, but not much larger)
1368* Reduce to bidiagonal form without LQ decomposition
1369*
1370 IE = 1
1371 ITAUQ = IE + M
1372 ITAUP = ITAUQ + M
1373 NWORK = ITAUP + M
1374*
1375* Bidiagonalize A
1376* Workspace: need 3*M [e, tauq, taup] + N [work]
1377* Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1378*
1379 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1380 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1381 $ IERR )
1382 IF( WNTQN ) THEN
1383*
1384* Path 5tn (N > M, JOBZ='N')
1385* Perform bidiagonal SVD, only computing singular values
1386* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1387*
1388 CALL DBDSDC( 'l', 'n', M, S, WORK( IE ), DUM, 1, DUM, 1,
1389 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1390 ELSE IF( WNTQO ) THEN
1391* Path 5to (N > M, JOBZ='O')
1392 LDWKVT = M
1393 IVT = NWORK
1394.GE. IF( LWORK M*N + 3*M + BDSPAC ) THEN
1395*
1396* WORK( IVT ) is M by N
1397*
1398 CALL DLASET( 'f', M, N, ZERO, ZERO, WORK( IVT ),
1399 $ LDWKVT )
1400 NWORK = IVT + LDWKVT*N
1401* IL is unused; silence compile warnings
1402 IL = -1
1403 ELSE
1404*
1405* WORK( IVT ) is M by M
1406*
1407 NWORK = IVT + LDWKVT*M
1408 IL = NWORK
1409*
1410* WORK(IL) is M by CHUNK
1411*
1412 CHUNK = ( LWORK - M*M - 3*M ) / M
1413 END IF
1414*
1415* Perform bidiagonal SVD, computing left singular vectors
1416* of bidiagonal matrix in U and computing right singular
1417* vectors of bidiagonal matrix in WORK(IVT)
1418* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1419*
1420 CALL DBDSDC( 'l', 'i', M, S, WORK( IE ), U, LDU,
1421 $ WORK( IVT ), LDWKVT, DUM, IDUM,
1422 $ WORK( NWORK ), IWORK, INFO )
1423*
1424* Overwrite U by left singular vectors of A
1425* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1426* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1427*
1428 CALL DORMBR( 'q', 'l', 'n', M, M, N, A, LDA,
1429 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1430 $ LWORK - NWORK + 1, IERR )
1431*
1432.GE. IF( LWORK M*N + 3*M + BDSPAC ) THEN
1433*
1434* Path 5to-fast
1435* Overwrite WORK(IVT) by left singular vectors of A
1436* Workspace: need 3*M [e, tauq, taup] + M*N [VT] + M [work]
1437* Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work]
1438*
1439 CALL DORMBR( 'p', 'r', 't', M, N, M, A, LDA,
1440 $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1441 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1442*
1443* Copy right singular vectors of A from WORK(IVT) to A
1444*
1445 CALL DLACPY( 'f', M, N, WORK( IVT ), LDWKVT, A, LDA )
1446 ELSE
1447*
1448* Path 5to-slow
1449* Generate P**T in A
1450* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work]
1451* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1452*
1453 CALL DORGBR( 'p', M, N, M, A, LDA, WORK( ITAUP ),
1454 $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1455*
1456* Multiply Q in A by right singular vectors of
1457* bidiagonal matrix in WORK(IVT), storing result in
1458* WORK(IL) and copying to A
1459* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M*NB [L]
1460* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N [L]
1461*
1462 DO 40 I = 1, N, CHUNK
1463 BLK = MIN( N - I + 1, CHUNK )
1464 CALL DGEMM( 'n', 'n', M, BLK, M, ONE, WORK( IVT ),
1465 $ LDWKVT, A( 1, I ), LDA, ZERO,
1466 $ WORK( IL ), M )
1467 CALL DLACPY( 'f', M, BLK, WORK( IL ), M, A( 1, I ),
1468 $ LDA )
1469 40 CONTINUE
1470 END IF
1471 ELSE IF( WNTQS ) THEN
1472*
1473* Path 5ts (N > M, JOBZ='S')
1474* Perform bidiagonal SVD, computing left singular vectors
1475* of bidiagonal matrix in U and computing right singular
1476* vectors of bidiagonal matrix in VT
1477* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1478*
1479 CALL DLASET( 'f', M, N, ZERO, ZERO, VT, LDVT )
1480 CALL DBDSDC( 'l', 'i', M, S, WORK( IE ), U, LDU, VT,
1481 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1482 $ INFO )
1483*
1484* Overwrite U by left singular vectors of A and VT
1485* by right singular vectors of A
1486* Workspace: need 3*M [e, tauq, taup] + M [work]
1487* Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1488*
1489 CALL DORMBR( 'q', 'l', 'n', M, M, N, A, LDA,
1490 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1491 $ LWORK - NWORK + 1, IERR )
1492 CALL DORMBR( 'p', 'r', 't', M, N, M, A, LDA,
1493 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1494 $ LWORK - NWORK + 1, IERR )
1495 ELSE IF( WNTQA ) THEN
1496*
1497* Path 5ta (N > M, JOBZ='A')
1498* Perform bidiagonal SVD, computing left singular vectors
1499* of bidiagonal matrix in U and computing right singular
1500* vectors of bidiagonal matrix in VT
1501* Workspace: need 3*M [e, tauq, taup] + BDSPAC
1502*
1503 CALL DLASET( 'f', N, N, ZERO, ZERO, VT, LDVT )
1504 CALL DBDSDC( 'l', 'i', M, S, WORK( IE ), U, LDU, VT,
1505 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1506 $ INFO )
1507*
1508* Set the right corner of VT to identity matrix
1509*
1510.GT. IF( NM ) THEN
1511 CALL DLASET( 'f', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1512 $ LDVT )
1513 END IF
1514*
1515* Overwrite U by left singular vectors of A and VT
1516* by right singular vectors of A
1517* Workspace: need 3*M [e, tauq, taup] + N [work]
1518* Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1519*
1520 CALL DORMBR( 'q', 'l', 'n', M, M, N, A, LDA,
1521 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1522 $ LWORK - NWORK + 1, IERR )
1523 CALL DORMBR( 'p', 'r', 't', N, N, M, A, LDA,
1524 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1525 $ LWORK - NWORK + 1, IERR )
1526 END IF
1527*
1528 END IF
1529*
1530 END IF
1531*
1532* Undo scaling if necessary
1533*
1534.EQ. IF( ISCL1 ) THEN
1535.GT. IF( ANRMBIGNUM )
1536 $ CALL DLASCL( 'g', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1537 $ IERR )
1538.LT. IF( ANRMSMLNUM )
1539 $ CALL DLASCL( 'g', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1540 $ IERR )
1541 END IF
1542*
1543* Return optimal workspace in WORK(1)
1544*
1545 WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
1546*
1547 RETURN
1548*
1549* End of DGESDD
1550*
1551 END
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dbdsdc
subroutine dbdsdc(uplo, compq, n, d, e, u, ldu, vt, ldvt, q, iq, work, iwork, info)
DBDSDC
Definition dbdsdc.f:205
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dorgbr
subroutine dorgbr(vect, m, n, k, a, lda, tau, work, lwork, info)
DORGBR
Definition dorgbr.f:157
dgebrd
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:205
dgelqf
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:143
dgeqrf
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
dgesdd
subroutine dgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
DGESDD
Definition dgesdd.f:219
dormbr
subroutine dormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMBR
Definition dormbr.f:195
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:127
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21