1 // Copyright ©2015 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "gonum.org/v1/gonum/blas"
11 "gonum.org/v1/gonum/blas/blas64"
12 "gonum.org/v1/gonum/lapack"
15 const noSVDO = "dgesvd: not coded for overwrite"
17 // Dgesvd computes the singular value decomposition of the input matrix A.
19 // The singular value decomposition is
20 // A = U * Sigma * V^T
21 // where Sigma is an m×n diagonal matrix containing the singular values of A,
22 // U is an m×m orthogonal matrix and V is an n×n orthogonal matrix. The first
23 // min(m,n) columns of U and V are the left and right singular vectors of A
26 // jobU and jobVT are options for computing the singular vectors. The behavior
28 // jobU == lapack.SVDAll All m columns of U are returned in u
29 // jobU == lapack.SVDInPlace The first min(m,n) columns are returned in u
30 // jobU == lapack.SVDOverwrite The first min(m,n) columns of U are written into a
31 // jobU == lapack.SVDNone The columns of U are not computed.
32 // The behavior is the same for jobVT and the rows of V^T. At most one of jobU
33 // and jobVT can equal lapack.SVDOverwrite, and Dgesvd will panic otherwise.
35 // On entry, a contains the data for the m×n matrix A. During the call to Dgesvd
36 // the data is overwritten. On exit, A contains the appropriate singular vectors
37 // if either job is lapack.SVDOverwrite.
39 // s is a slice of length at least min(m,n) and on exit contains the singular
40 // values in decreasing order.
42 // u contains the left singular vectors on exit, stored column-wise. If
43 // jobU == lapack.SVDAll, u is of size m×m. If jobU == lapack.SVDInPlace u is
44 // of size m×min(m,n). If jobU == lapack.SVDOverwrite or lapack.SVDNone, u is
47 // vt contains the left singular vectors on exit, stored row-wise. If
48 // jobV == lapack.SVDAll, vt is of size n×m. If jobVT == lapack.SVDInPlace vt is
49 // of size min(m,n)×n. If jobVT == lapack.SVDOverwrite or lapack.SVDNone, vt is
52 // work is a slice for storing temporary memory, and lwork is the usable size of
53 // the slice. lwork must be at least max(5*min(m,n), 3*min(m,n)+max(m,n)).
54 // If lwork == -1, instead of performing Dgesvd, the optimal work length will be
55 // stored into work[0]. Dgesvd will panic if the working memory has insufficient
58 // Dgesvd returns whether the decomposition successfully completed.
59 func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool) {
61 checkMatrix(m, n, a, lda)
62 if jobU == lapack.SVDAll {
63 checkMatrix(m, m, u, ldu)
64 } else if jobU == lapack.SVDInPlace {
65 checkMatrix(m, minmn, u, ldu)
67 if jobVT == lapack.SVDAll {
68 checkMatrix(n, n, vt, ldvt)
69 } else if jobVT == lapack.SVDInPlace {
70 checkMatrix(minmn, n, vt, ldvt)
72 if jobU == lapack.SVDOverwrite && jobVT == lapack.SVDOverwrite {
73 panic("lapack: both jobU and jobVT are lapack.SVDOverwrite")
78 if jobU == lapack.SVDOverwrite || jobVT == lapack.SVDOverwrite {
85 wantua := jobU == lapack.SVDAll
86 wantus := jobU == lapack.SVDInPlace
87 wantuas := wantua || wantus
88 wantuo := jobU == lapack.SVDOverwrite
89 wantun := jobU == lapack.None
91 wantva := jobVT == lapack.SVDAll
92 wantvs := jobVT == lapack.SVDInPlace
93 wantvas := wantva || wantvs
94 wantvo := jobVT == lapack.SVDOverwrite
95 wantvn := jobVT == lapack.None
97 bi := blas64.Implementation()
100 // Compute optimal space for subroutines.
102 opts := string(jobU) + string(jobVT)
103 var wrkbl, bdspac int
105 mnthr = impl.Ilaenv(6, "DGESVD", opts, m, n, 0, 0)
107 impl.Dgeqrf(m, n, a, lda, nil, work, -1)
108 lwork_dgeqrf := int(work[0])
109 impl.Dorgqr(m, n, n, a, lda, nil, work, -1)
110 lwork_dorgqr_n := int(work[0])
111 impl.Dorgqr(m, m, n, a, lda, nil, work, -1)
112 lwork_dorgqr_m := int(work[0])
113 impl.Dgebrd(n, n, a, lda, s, nil, nil, nil, work, -1)
114 lwork_dgebrd := int(work[0])
115 impl.Dorgbr(lapack.ApplyP, n, n, n, a, lda, nil, work, -1)
116 lwork_dorgbr_p := int(work[0])
117 impl.Dorgbr(lapack.ApplyQ, n, n, n, a, lda, nil, work, -1)
118 lwork_dorgbr_q := int(work[0])
124 maxwrk = n + lwork_dgeqrf
125 maxwrk = max(maxwrk, 3*n+lwork_dgebrd)
126 if wantvo || wantvas {
127 maxwrk = max(maxwrk, 3*n+lwork_dorgbr_p)
129 maxwrk = max(maxwrk, bdspac)
130 } else if wantuo && wantvn {
132 wrkbl = n + lwork_dgeqrf
133 wrkbl = max(wrkbl, n+lwork_dorgqr_n)
134 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
135 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
136 wrkbl = max(wrkbl, bdspac)
137 maxwrk = max(n*n+wrkbl, n*n+m*n+n)
138 } else if wantuo && wantvs {
140 wrkbl = n + lwork_dgeqrf
141 wrkbl = max(wrkbl, n+lwork_dorgqr_n)
142 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
143 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
144 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
145 wrkbl = max(wrkbl, bdspac)
146 maxwrk = max(n*n+wrkbl, n*n+m*n+n)
147 } else if wantus && wantvn {
149 wrkbl = n + lwork_dgeqrf
150 wrkbl = max(wrkbl, n+lwork_dorgqr_n)
151 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
152 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
153 wrkbl = max(wrkbl, bdspac)
155 } else if wantus && wantvo {
157 wrkbl = n + lwork_dgeqrf
158 wrkbl = max(wrkbl, n+lwork_dorgqr_n)
159 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
160 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
161 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
162 wrkbl = max(wrkbl, bdspac)
163 maxwrk = 2*n*n + wrkbl
164 } else if wantus && wantvas {
166 wrkbl = n + lwork_dgeqrf
167 wrkbl = max(wrkbl, n+lwork_dorgqr_n)
168 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
169 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
170 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
171 wrkbl = max(wrkbl, bdspac)
173 } else if wantua && wantvn {
175 wrkbl = n + lwork_dgeqrf
176 wrkbl = max(wrkbl, n+lwork_dorgqr_m)
177 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
178 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
179 wrkbl = max(wrkbl, bdspac)
181 } else if wantua && wantvo {
183 wrkbl = n + lwork_dgeqrf
184 wrkbl = max(wrkbl, n+lwork_dorgqr_m)
185 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
186 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
187 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
188 wrkbl = max(wrkbl, bdspac)
189 maxwrk = 2*n*n + wrkbl
190 } else if wantua && wantvas {
192 wrkbl = n + lwork_dgeqrf
193 wrkbl = max(wrkbl, n+lwork_dorgqr_m)
194 wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
195 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
196 wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
197 wrkbl = max(wrkbl, bdspac)
202 impl.Dgebrd(m, n, a, lda, s, nil, nil, nil, work, -1)
203 lwork_dgebrd := int(work[0])
204 maxwrk = 3*n + lwork_dgebrd
205 if wantus || wantuo {
206 impl.Dorgbr(lapack.ApplyQ, m, n, n, a, lda, nil, work, -1)
207 lwork_dorgbr_q = int(work[0])
208 maxwrk = max(maxwrk, 3*n+lwork_dorgbr_q)
211 impl.Dorgbr(lapack.ApplyQ, m, m, n, a, lda, nil, work, -1)
212 lwork_dorgbr_q := int(work[0])
213 maxwrk = max(maxwrk, 3*n+lwork_dorgbr_q)
216 maxwrk = max(maxwrk, 3*n+lwork_dorgbr_p)
218 maxwrk = max(maxwrk, bdspac)
221 mnthr = impl.Ilaenv(6, "DGESVD", opts, m, n, 0, 0)
224 impl.Dgelqf(m, n, a, lda, nil, work, -1)
225 lwork_dgelqf := int(work[0])
226 impl.Dorglq(n, n, m, nil, n, nil, work, -1)
227 lwork_dorglq_n := int(work[0])
228 impl.Dorglq(m, n, m, a, lda, nil, work, -1)
229 lwork_dorglq_m := int(work[0])
230 impl.Dgebrd(m, m, a, lda, s, nil, nil, nil, work, -1)
231 lwork_dgebrd := int(work[0])
232 impl.Dorgbr(lapack.ApplyP, m, m, m, a, n, nil, work, -1)
233 lwork_dorgbr_p := int(work[0])
234 impl.Dorgbr(lapack.ApplyQ, m, m, m, a, n, nil, work, -1)
235 lwork_dorgbr_q := int(work[0])
240 maxwrk = m + lwork_dgelqf
241 maxwrk = max(maxwrk, 3*m+lwork_dgebrd)
242 if wantuo || wantuas {
243 maxwrk = max(maxwrk, 3*m+lwork_dorgbr_q)
245 maxwrk = max(maxwrk, bdspac)
246 } else if wantvo && wantun {
248 wrkbl = m + lwork_dgelqf
249 wrkbl = max(wrkbl, m+lwork_dorglq_m)
250 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
251 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
252 wrkbl = max(wrkbl, bdspac)
253 maxwrk = max(m*m+wrkbl, m*m+m*n+m)
254 } else if wantvo && wantuas {
256 wrkbl = m + lwork_dgelqf
257 wrkbl = max(wrkbl, m+lwork_dorglq_m)
258 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
259 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
260 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
261 wrkbl = max(wrkbl, bdspac)
262 maxwrk = max(m*m+wrkbl, m*m+m*n+m)
263 } else if wantvs && wantun {
265 wrkbl = m + lwork_dgelqf
266 wrkbl = max(wrkbl, m+lwork_dorglq_m)
267 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
268 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
269 wrkbl = max(wrkbl, bdspac)
271 } else if wantvs && wantuo {
273 wrkbl = m + lwork_dgelqf
274 wrkbl = max(wrkbl, m+lwork_dorglq_m)
275 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
276 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
277 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
278 wrkbl = max(wrkbl, bdspac)
279 maxwrk = 2*m*m + wrkbl
280 } else if wantvs && wantuas {
282 wrkbl = m + lwork_dgelqf
283 wrkbl = max(wrkbl, m+lwork_dorglq_m)
284 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
285 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
286 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
287 wrkbl = max(wrkbl, bdspac)
289 } else if wantva && wantun {
291 wrkbl = m + lwork_dgelqf
292 wrkbl = max(wrkbl, m+lwork_dorglq_n)
293 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
294 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
295 wrkbl = max(wrkbl, bdspac)
297 } else if wantva && wantuo {
299 wrkbl = m + lwork_dgelqf
300 wrkbl = max(wrkbl, m+lwork_dorglq_n)
301 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
302 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
303 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
304 wrkbl = max(wrkbl, bdspac)
305 maxwrk = 2*m*m + wrkbl
306 } else if wantva && wantuas {
308 wrkbl = m + lwork_dgelqf
309 wrkbl = max(wrkbl, m+lwork_dorglq_n)
310 wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
311 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
312 wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
313 wrkbl = max(wrkbl, bdspac)
318 impl.Dgebrd(m, n, a, lda, s, nil, nil, nil, work, -1)
319 lwork_dgebrd = int(work[0])
320 maxwrk = 3*m + lwork_dgebrd
321 if wantvs || wantvo {
322 impl.Dorgbr(lapack.ApplyP, m, n, m, a, n, nil, work, -1)
323 lwork_dorgbr_p = int(work[0])
324 maxwrk = max(maxwrk, 3*m+lwork_dorgbr_p)
327 impl.Dorgbr(lapack.ApplyP, n, n, m, a, n, nil, work, -1)
328 lwork_dorgbr_p = int(work[0])
329 maxwrk = max(maxwrk, 3*m+lwork_dorgbr_p)
332 maxwrk = max(maxwrk, 3*m+lwork_dorgbr_q)
334 maxwrk = max(maxwrk, bdspac)
338 minWork := max(1, 5*minmn)
339 if !((wantun && m >= mnthr) || (wantvn && n >= mnthr)) {
340 minWork = max(minWork, 3*minmn+max(m, n))
344 if len(work) < lwork {
351 if m == 0 || n == 0 {
355 maxwrk = max(maxwrk, minWork)
356 work[0] = float64(maxwrk)
361 // Perform decomposition.
363 smlnum := math.Sqrt(dlamchS) / eps
366 // Scale A if max element outside range [smlnum, bignum].
367 anrm := impl.Dlange(lapack.MaxAbs, m, n, a, lda, nil)
369 if anrm > 0 && anrm < smlnum {
371 impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, m, n, a, lda)
372 } else if anrm > bignum {
374 impl.Dlascl(lapack.General, 0, 0, anrm, bignum, m, n, a, lda)
379 // If A has sufficiently more rows than columns, use the QR decomposition.
387 // Compute A = Q * R.
388 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
391 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
396 // Bidiagonalize R in A.
397 impl.Dgebrd(n, n, a, lda, s, work[ie:], work[itauq:],
398 work[itaup:], work[iwork:], lwork-iwork)
400 if wantvo || wantvas {
402 impl.Dorgbr(lapack.ApplyP, n, n, n, a, lda, work[itaup:],
403 work[iwork:], lwork-iwork)
408 // Perform bidiagonal QR iteration computing right singular vectors
409 // of A in A if desired.
410 ok = impl.Dbdsqr(blas.Upper, n, ncvt, 0, 0, s, work[ie:],
411 a, lda, work, 1, work, 1, work[iwork:])
413 // If right singular vectors desired in VT, copy them there.
415 impl.Dlacpy(blas.All, n, n, a, lda, vt, ldvt)
417 } else if wantuo && wantvn {
420 } else if wantuo && wantvas {
426 if lwork >= n*n+max(4*n, bdspac) {
427 // Sufficient workspace for a fast algorithm.
430 if lwork >= wrkbl+lda*n {
435 itau := ir + ldworkr*n
437 // Compute A = Q * R.
438 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
440 // Copy R to work[ir:], zeroing out below it.
441 impl.Dlacpy(blas.Upper, n, n, a, lda, work[ir:], ldworkr)
442 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[ir+ldworkr:], ldworkr)
445 impl.Dorgqr(m, n, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
451 // Bidiagonalize R in work[ir:].
452 impl.Dgebrd(n, n, work[ir:], ldworkr, s, work[ie:],
453 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
455 // Generate left vectors bidiagonalizing R in work[ir:].
456 impl.Dorgbr(lapack.ApplyQ, n, n, n, work[ir:], ldworkr,
457 work[itauq:], work[iwork:], lwork-iwork)
460 // Perform bidiagonal QR iteration, compuing left singular
461 // vectors of R in work[ir:].
462 ok = impl.Dbdsqr(blas.Upper, n, 0, n, 0, s, work[ie:], work, 1,
463 work[ir:], ldworkr, work, 1, work[iwork:])
465 // Multiply Q in A by left singular vectors of R in
466 // work[ir:], storing result in U.
467 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, a, lda,
468 work[ir:], ldworkr, 0, u, ldu)
470 // Insufficient workspace for a fast algorithm.
474 // Compute A = Q*R, copying result to U.
475 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
476 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
479 impl.Dorgqr(m, n, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
485 // Zero out below R in A.
486 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
488 // Bidiagonalize R in A.
489 impl.Dgebrd(n, n, a, lda, s, work[ie:],
490 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
492 // Multiply Q in U by left vectors bidiagonalizing R.
493 impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
494 a, lda, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
497 // Perform bidiagonal QR iteration, computing left
498 // singular vectors of A in U.
499 ok = impl.Dbdsqr(blas.Upper, n, 0, m, 0, s, work[ie:], work, 1,
500 u, ldu, work, 1, work[iwork:])
507 if lwork >= n*n+max(4*n, bdspac) {
508 // Sufficient workspace for a fast algorithm.
511 if lwork >= wrkbl+lda*n {
516 itau := iu + ldworku*n
519 // Compute A = Q * R.
520 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
521 // Copy R to work[iu:], zeroing out below it.
522 impl.Dlacpy(blas.Upper, n, n, a, lda, work[iu:], ldworku)
523 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[iu+ldworku:], ldworku)
526 impl.Dorgqr(m, n, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
533 // Bidiagonalize R in work[iu:], copying result to VT.
534 impl.Dgebrd(n, n, work[iu:], ldworku, s, work[ie:],
535 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
536 impl.Dlacpy(blas.Upper, n, n, work[iu:], ldworku, vt, ldvt)
538 // Generate left bidiagonalizing vectors in work[iu:].
539 impl.Dorgbr(lapack.ApplyQ, n, n, n, work[iu:], ldworku,
540 work[itauq:], work[iwork:], lwork-iwork)
542 // Generate right bidiagonalizing vectors in VT.
543 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
544 work[itaup:], work[iwork:], lwork-iwork)
547 // Perform bidiagonal QR iteration, computing left singular
548 // vectors of R in work[iu:], and computing right singular
549 // vectors of R in VT.
550 ok = impl.Dbdsqr(blas.Upper, n, n, n, 0, s, work[ie:],
551 vt, ldvt, work[iu:], ldworku, work, 1, work[iwork:])
553 // Multiply Q in A by left singular vectors of R in
554 // work[iu:], storing result in U.
555 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, a, lda,
556 work[iu:], ldworku, 0, u, ldu)
558 // Insufficient workspace for a fast algorithm.
562 // Compute A = Q * R, copying result to U.
563 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
564 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
567 impl.Dorgqr(m, n, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
569 // Copy R to VT, zeroing out below it.
570 impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
571 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, vt[ldvt:], ldvt)
578 // Bidiagonalize R in VT.
579 impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
580 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
582 // Multiply Q in U by left bidiagonalizing vectors in VT.
583 impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
584 vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
586 // Generate right bidiagonalizing vectors in VT.
587 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
588 work[itaup:], work[iwork:], lwork-iwork)
591 // Perform bidiagonal QR iteration, computing left singular
592 // vectors of A in U and computing right singular vectors
594 ok = impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
595 vt, ldvt, u, ldu, work, 1, work[iwork:])
601 if lwork >= n*n+max(max(n+m, 4*n), bdspac) {
602 // Sufficient workspace for a fast algorithm.
605 if lwork >= wrkbl+lda*n {
610 itau := ir + ldworkr*n
613 // Compute A = Q*R, copying result to U.
614 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
615 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
617 // Copy R to work[ir:], zeroing out below it.
618 impl.Dlacpy(blas.Upper, n, n, a, lda, work[ir:], ldworkr)
619 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[ir+ldworkr:], ldworkr)
622 impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
628 // Bidiagonalize R in work[ir:].
629 impl.Dgebrd(n, n, work[ir:], ldworkr, s, work[ie:],
630 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
632 // Generate left bidiagonalizing vectors in work[ir:].
633 impl.Dorgbr(lapack.ApplyQ, n, n, n, work[ir:], ldworkr,
634 work[itauq:], work[iwork:], lwork-iwork)
637 // Perform bidiagonal QR iteration, computing left singular
638 // vectors of R in work[ir:].
639 ok = impl.Dbdsqr(blas.Upper, n, 0, n, 0, s, work[ie:], work, 1,
640 work[ir:], ldworkr, work, 1, work[iwork:])
642 // Multiply Q in U by left singular vectors of R in
643 // work[ir:], storing result in A.
644 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, u, ldu,
645 work[ir:], ldworkr, 0, a, lda)
647 // Copy left singular vectors of A from A to U.
648 impl.Dlacpy(blas.All, m, n, a, lda, u, ldu)
650 // Insufficient workspace for a fast algorithm.
654 // Compute A = Q*R, copying result to U.
655 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
656 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
659 impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
665 // Zero out below R in A.
666 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
668 // Bidiagonalize R in A.
669 impl.Dgebrd(n, n, a, lda, s, work[ie:],
670 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
672 // Multiply Q in U by left bidiagonalizing vectors in A.
673 impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
674 a, lda, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
677 // Perform bidiagonal QR iteration, computing left
678 // singular vectors of A in U.
679 ok = impl.Dbdsqr(blas.Upper, n, 0, m, 0, s, work[ie:],
680 work, 1, u, ldu, work, 1, work[iwork:])
687 if lwork >= n*n+max(max(n+m, 4*n), bdspac) {
688 // Sufficient workspace for a fast algorithm.
691 if lwork >= wrkbl+lda*n {
696 itau := iu + ldworku*n
699 // Compute A = Q * R, copying result to U.
700 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
701 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
704 impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
706 // Copy R to work[iu:], zeroing out below it.
707 impl.Dlacpy(blas.Upper, n, n, a, lda, work[iu:], ldworku)
708 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[iu+ldworku:], ldworku)
715 // Bidiagonalize R in work[iu:], copying result to VT.
716 impl.Dgebrd(n, n, work[iu:], ldworku, s, work[ie:],
717 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
718 impl.Dlacpy(blas.Upper, n, n, work[iu:], ldworku, vt, ldvt)
720 // Generate left bidiagonalizing vectors in work[iu:].
721 impl.Dorgbr(lapack.ApplyQ, n, n, n, work[iu:], ldworku,
722 work[itauq:], work[iwork:], lwork-iwork)
724 // Generate right bidiagonalizing vectors in VT.
725 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
726 work[itaup:], work[iwork:], lwork-iwork)
729 // Perform bidiagonal QR iteration, computing left singular
730 // vectors of R in work[iu:] and computing right
731 // singular vectors of R in VT.
732 ok = impl.Dbdsqr(blas.Upper, n, n, n, 0, s, work[ie:],
733 vt, ldvt, work[iu:], ldworku, work, 1, work[iwork:])
735 // Multiply Q in U by left singular vectors of R in
736 // work[iu:], storing result in A.
737 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1,
738 u, ldu, work[iu:], ldworku, 0, a, lda)
740 // Copy left singular vectors of A from A to U.
741 impl.Dlacpy(blas.All, m, n, a, lda, u, ldu)
744 // Bidiagonalize R in VT.
745 impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
746 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
748 // Multiply Q in U by left bidiagonalizing vectors in VT.
749 impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans,
750 m, n, n, vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
752 // Generate right bidiagonalizing vectors in VT.
753 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
754 work[itaup:], work[iwork:], lwork-iwork)
757 // Perform bidiagonal QR iteration, computing left singular
758 // vectors of A in U and computing right singular vectors
760 ok = impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
761 vt, ldvt, u, ldu, work, 1, work[iwork:])
764 // Insufficient workspace for a fast algorithm.
768 // Compute A = Q*R, copying result to U.
769 impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
770 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
773 impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
775 // Copy R from A to VT, zeroing out below it.
776 impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
777 impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, vt[ldvt:], ldvt)
784 // Bidiagonalize R in VT.
785 impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
786 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
788 // Multiply Q in U by left bidiagonalizing vectors in VT.
789 impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans,
790 m, n, n, vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
792 // Generate right bidiagonizing vectors in VT.
793 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
794 work[itaup:], work[iwork:], lwork-iwork)
797 // Perform bidiagonal QR iteration, computing left singular
798 // vectors of A in U and computing right singular vectors
800 impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
801 vt, ldvt, u, ldu, work, 1, work[iwork:])
807 // M at least N, but not much larger.
814 impl.Dgebrd(m, n, a, lda, s, work[ie:], work[itauq:],
815 work[itaup:], work[iwork:], lwork-iwork)
817 // Left singular vectors are desired in U. Copy result to U and
818 // generate left biadiagonalizing vectors in U.
819 impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
827 impl.Dorgbr(lapack.ApplyQ, m, ncu, n, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
830 // Right singular vectors are desired in VT. Copy result to VT and
831 // generate left biadiagonalizing vectors in VT.
832 impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
833 impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt, work[itaup:], work[iwork:], lwork-iwork)
843 if wantuas || wantuo {
849 if wantvas || wantvo {
855 if !wantuo && !wantvo {
856 // Perform bidiagonal QR iteration, if desired, computing left
857 // singular vectors in U and right singular vectors in VT.
858 ok = impl.Dbdsqr(blas.Upper, n, ncvt, nru, 0, s, work[ie:],
859 vt, ldvt, u, ldu, work, 1, work[iwork:])
861 // There will be two branches when the implementation is complete.
866 // A has more columns than rows. If A has sufficiently more columns than
867 // rows, first reduce using the LQ decomposition.
876 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
879 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
885 // Bidiagonalize L in A.
886 impl.Dgebrd(m, m, a, lda, s, work[ie:itauq],
887 work[itauq:itaup], work[itaup:iwork], work[iwork:], lwork-iwork)
888 if wantuo || wantuas {
889 impl.Dorgbr(lapack.ApplyQ, m, m, m, a, lda,
890 work[itauq:], work[iwork:], lwork-iwork)
894 if wantuo || wantuas {
898 // Perform bidiagonal QR iteration, computing left singular vectors
899 // of A in A if desired.
900 ok = impl.Dbdsqr(blas.Upper, m, 0, nru, 0, s, work[ie:],
901 work, 1, a, lda, work, 1, work[iwork:])
903 // If left singular vectors desired in U, copy them there.
905 impl.Dlacpy(blas.All, m, m, a, lda, u, ldu)
907 } else if wantvo && wantun {
910 } else if wantvo && wantuas {
916 if lwork >= m*m+max(4*m, bdspac) {
917 // Sufficient workspace for a fast algorithm.
920 if lwork >= wrkbl+lda*m {
925 itau := ir + ldworkr*m
929 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
931 // Copy L to work[ir:], zeroing out above it.
932 impl.Dlacpy(blas.Lower, m, m, a, lda, work[ir:], ldworkr)
933 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[ir+1:], ldworkr)
936 impl.Dorglq(m, n, m, a, lda, work[itau:], work[iwork:], lwork-iwork)
942 // Bidiagonalize L in work[ir:].
943 impl.Dgebrd(m, m, work[ir:], ldworkr, s, work[ie:],
944 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
946 // Generate right vectors bidiagonalizing L in work[ir:].
947 impl.Dorgbr(lapack.ApplyP, m, m, m, work[ir:], ldworkr,
948 work[itaup:], work[iwork:], lwork-iwork)
951 // Perform bidiagonal QR iteration, computing right singular
952 // vectors of L in work[ir:].
953 ok = impl.Dbdsqr(blas.Upper, m, m, 0, 0, s, work[ie:],
954 work[ir:], ldworkr, work, 1, work, 1, work[iwork:])
956 // Multiply right singular vectors of L in work[ir:] by
957 // Q in A, storing result in VT.
958 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
959 work[ir:], ldworkr, a, lda, 0, vt, ldvt)
961 // Insufficient workspace for a fast algorithm.
966 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
968 // Copy result to VT.
969 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
972 impl.Dorglq(m, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
978 // Zero out above L in A.
979 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
981 // Bidiagonalize L in A.
982 impl.Dgebrd(m, m, a, lda, s, work[ie:],
983 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
985 // Multiply right vectors bidiagonalizing L by Q in VT.
986 impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
987 a, lda, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
990 // Perform bidiagonal QR iteration, computing right
991 // singular vectors of A in VT.
992 ok = impl.Dbdsqr(blas.Upper, m, n, 0, 0, s, work[ie:],
993 vt, ldvt, work, 1, work, 1, work[iwork:])
1000 if lwork >= m*m+max(4*m, bdspac) {
1001 // Sufficient workspace for a fast algorithm.
1004 if lwork >= wrkbl+lda*m {
1009 itau := iu + ldworku*m
1013 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1015 // Copy L to work[iu:], zeroing out above it.
1016 impl.Dlacpy(blas.Lower, m, m, a, lda, work[iu:], ldworku)
1017 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[iu+1:], ldworku)
1020 impl.Dorglq(m, n, m, a, lda, work[itau:], work[iwork:], lwork-iwork)
1026 // Bidiagonalize L in work[iu:], copying result to U.
1027 impl.Dgebrd(m, m, work[iu:], ldworku, s, work[ie:],
1028 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
1029 impl.Dlacpy(blas.Lower, m, m, work[iu:], ldworku, u, ldu)
1031 // Generate right bidiagionalizing vectors in work[iu:].
1032 impl.Dorgbr(lapack.ApplyP, m, m, m, work[iu:], ldworku,
1033 work[itaup:], work[iwork:], lwork-iwork)
1035 // Generate left bidiagonalizing vectors in U.
1036 impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
1039 // Perform bidiagonal QR iteration, computing left singular
1040 // vectors of L in U and computing right singular vectors of
1042 ok = impl.Dbdsqr(blas.Upper, m, m, m, 0, s, work[ie:],
1043 work[iu:], ldworku, u, ldu, work, 1, work[iwork:])
1045 // Multiply right singular vectors of L in work[iu:] by
1046 // Q in A, storing result in VT.
1047 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
1048 work[iu:], ldworku, a, lda, 0, vt, ldvt)
1050 // Insufficient workspace for a fast algorithm.
1054 // Compute A = L*Q, copying result to VT.
1055 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1056 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1058 // Generate Q in VT.
1059 impl.Dorglq(m, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
1061 // Copy L to U, zeroing out above it.
1062 impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
1063 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, u[1:], ldu)
1070 // Bidiagonalize L in U.
1071 impl.Dgebrd(m, m, u, ldu, s, work[ie:],
1072 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
1074 // Multiply right bidiagonalizing vectors in U by Q in VT.
1075 impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
1076 u, ldu, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
1078 // Generate left bidiagonalizing vectors in U.
1079 impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
1082 // Perform bidiagonal QR iteration, computing left singular
1083 // vectors of A in U and computing right singular vectors
1085 impl.Dbdsqr(blas.Upper, m, n, m, 0, s, work[ie:], vt, ldvt,
1086 u, ldu, work, 1, work[iwork:])
1092 if lwork >= m*m+max(max(n+m, 4*m), bdspac) {
1093 // Sufficient workspace for a fast algorithm.
1096 if lwork >= wrkbl+lda*m {
1101 itau := ir + ldworkr*m
1104 // Compute A = L*Q, copying result to VT.
1105 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1106 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1108 // Copy L to work[ir:], zeroing out above it.
1109 impl.Dlacpy(blas.Lower, m, m, a, lda, work[ir:], ldworkr)
1110 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[ir+1:], ldworkr)
1112 // Generate Q in VT.
1113 impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
1120 // Bidiagonalize L in work[ir:].
1121 impl.Dgebrd(m, m, work[ir:], ldworkr, s, work[ie:],
1122 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
1124 // Generate right bidiagonalizing vectors in work[ir:].
1125 impl.Dorgbr(lapack.ApplyP, m, m, m, work[ir:], ldworkr,
1126 work[itaup:], work[iwork:], lwork-iwork)
1129 // Perform bidiagonal QR iteration, computing right
1130 // singular vectors of L in work[ir:].
1131 ok = impl.Dbdsqr(blas.Upper, m, m, 0, 0, s, work[ie:],
1132 work[ir:], ldworkr, work, 1, work, 1, work[iwork:])
1134 // Multiply right singular vectors of L in work[ir:] by
1135 // Q in VT, storing result in A.
1136 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
1137 work[ir:], ldworkr, vt, ldvt, 0, a, lda)
1139 // Copy right singular vectors of A from A to VT.
1140 impl.Dlacpy(blas.All, m, n, a, lda, vt, ldvt)
1142 // Insufficient workspace for a fast algorithm.
1145 // Compute A = L * Q, copying result to VT.
1146 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1147 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1149 // Generate Q in VT.
1150 impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
1157 // Zero out above L in A.
1158 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
1160 // Bidiagonalize L in A.
1161 impl.Dgebrd(m, m, a, lda, s, work[ie:], work[itauq:],
1162 work[itaup:], work[iwork:], lwork-iwork)
1164 // Multiply right bidiagonalizing vectors in A by Q in VT.
1165 impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
1166 a, lda, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
1169 // Perform bidiagonal QR iteration, computing right singular
1170 // vectors of A in VT.
1171 ok = impl.Dbdsqr(blas.Upper, m, n, 0, 0, s, work[ie:],
1172 vt, ldvt, work, 1, work, 1, work[iwork:])
1178 if lwork >= m*m+max(max(m+n, 4*m), bdspac) {
1179 // Sufficient workspace for a fast algorithm.
1183 if lwork >= wrkbl+lda*m {
1188 itau := iu + ldworku*m
1191 // Generate A = L * Q copying result to VT.
1192 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1193 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1195 // Generate Q in VT.
1196 impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
1198 // Copy L to work[iu:], zeroing out above it.
1199 impl.Dlacpy(blas.Lower, m, m, a, lda, work[iu:], ldworku)
1200 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[iu+1:], ldworku)
1206 // Bidiagonalize L in work[iu:], copying result to U.
1207 impl.Dgebrd(m, m, work[iu:], ldworku, s, work[ie:],
1208 work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
1209 impl.Dlacpy(blas.Lower, m, m, work[iu:], ldworku, u, ldu)
1211 // Generate right bidiagonalizing vectors in work[iu:].
1212 impl.Dorgbr(lapack.ApplyP, m, m, m, work[iu:], ldworku,
1213 work[itaup:], work[iwork:], lwork-iwork)
1215 // Generate left bidiagonalizing vectors in U.
1216 impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
1219 // Perform bidiagonal QR iteration, computing left singular
1220 // vectors of L in U and computing right singular vectors
1221 // of L in work[iu:].
1222 ok = impl.Dbdsqr(blas.Upper, m, m, m, 0, s, work[ie:],
1223 work[iu:], ldworku, u, ldu, work, 1, work[iwork:])
1225 // Multiply right singular vectors of L in work[iu:]
1226 // Q in VT, storing result in A.
1227 bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
1228 work[iu:], ldworku, vt, ldvt, 0, a, lda)
1230 // Copy right singular vectors of A from A to VT.
1231 impl.Dlacpy(blas.All, m, n, a, lda, vt, ldvt)
1233 // Insufficient workspace for a fast algorithm.
1237 // Compute A = L * Q, copying result to VT.
1238 impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
1239 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1241 // Generate Q in VT.
1242 impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
1244 // Copy L to U, zeroing out above it.
1245 impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
1246 impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, u[1:], ldu)
1253 // Bidiagonalize L in U.
1254 impl.Dgebrd(m, m, u, ldu, s, work[ie:], work[itauq:],
1255 work[itaup:], work[iwork:], lwork-iwork)
1257 // Multiply right bidiagonalizing vectors in U by Q in VT.
1258 impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
1259 u, ldu, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
1261 // Generate left bidiagonalizing vectors in U.
1262 impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
1265 // Perform bidiagonal QR iteration, computing left singular
1266 // vectors of A in U and computing right singular vectors
1268 ok = impl.Dbdsqr(blas.Upper, m, n, m, 0, s, work[ie:],
1269 vt, ldvt, u, ldu, work, 1, work[iwork:])
1275 // N at least M, but not much larger.
1282 impl.Dgebrd(m, n, a, lda, s, work[ie:], work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
1284 // If left singular vectors desired in U, copy result to U and
1285 // generate left bidiagonalizing vectors in U.
1286 impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
1287 impl.Dorgbr(lapack.ApplyQ, m, m, n, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
1290 // If right singular vectors desired in VT, copy result to VT
1291 // and generate right bidiagonalizing vectors in VT.
1292 impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
1299 impl.Dorgbr(lapack.ApplyP, nrvt, n, m, vt, ldvt, work[itaup:], work[iwork:], lwork-iwork)
1309 if wantuas || wantuo {
1312 if wantvas || wantvo {
1315 if !wantuo && !wantvo {
1316 // Perform bidiagonal QR iteration, if desired, computing left
1317 // singular vectors in U and computing right singular vectors in
1319 ok = impl.Dbdsqr(blas.Lower, m, ncvt, nru, 0, s, work[ie:],
1320 vt, ldvt, u, ldu, work, 1, work[iwork:])
1322 // There will be two branches when the implementation is complete.
1329 for i := 0; i < minmn-1; i++ {
1330 work[i+1] = work[i+ie]
1334 for i := minmn - 2; i >= 0; i-- {
1335 work[i+1] = work[i+ie]
1339 // Undo scaling if necessary.
1342 impl.Dlascl(lapack.General, 0, 0, bignum, anrm, minmn, 1, s, minmn)
1344 if !ok && anrm > bignum {
1345 impl.Dlascl(lapack.General, 0, 0, bignum, anrm, minmn-1, 1, work[minmn:], minmn)
1348 impl.Dlascl(lapack.General, 0, 0, smlnum, anrm, minmn, 1, s, minmn)
1350 if !ok && anrm < smlnum {
1351 impl.Dlascl(lapack.General, 0, 0, smlnum, anrm, minmn-1, 1, work[minmn:], minmn)
1354 work[0] = float64(maxwrk)