--- /dev/null
+// Copyright ©2015 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package gonum
+
+import (
+ "math"
+
+ "gonum.org/v1/gonum/blas"
+ "gonum.org/v1/gonum/blas/blas64"
+ "gonum.org/v1/gonum/lapack"
+)
+
+const noSVDO = "dgesvd: not coded for overwrite"
+
+// Dgesvd computes the singular value decomposition of the input matrix A.
+//
+// The singular value decomposition is
+// A = U * Sigma * V^T
+// where Sigma is an m×n diagonal matrix containing the singular values of A,
+// U is an m×m orthogonal matrix and V is an n×n orthogonal matrix. The first
+// min(m,n) columns of U and V are the left and right singular vectors of A
+// respectively.
+//
+// jobU and jobVT are options for computing the singular vectors. The behavior
+// is as follows
+// jobU == lapack.SVDAll All m columns of U are returned in u
+// jobU == lapack.SVDInPlace The first min(m,n) columns are returned in u
+// jobU == lapack.SVDOverwrite The first min(m,n) columns of U are written into a
+// jobU == lapack.SVDNone The columns of U are not computed.
+// The behavior is the same for jobVT and the rows of V^T. At most one of jobU
+// and jobVT can equal lapack.SVDOverwrite, and Dgesvd will panic otherwise.
+//
+// On entry, a contains the data for the m×n matrix A. During the call to Dgesvd
+// the data is overwritten. On exit, A contains the appropriate singular vectors
+// if either job is lapack.SVDOverwrite.
+//
+// s is a slice of length at least min(m,n) and on exit contains the singular
+// values in decreasing order.
+//
+// u contains the left singular vectors on exit, stored column-wise. If
+// jobU == lapack.SVDAll, u is of size m×m. If jobU == lapack.SVDInPlace u is
+// of size m×min(m,n). If jobU == lapack.SVDOverwrite or lapack.SVDNone, u is
+// not used.
+//
+// vt contains the left singular vectors on exit, stored row-wise. If
+// jobV == lapack.SVDAll, vt is of size n×m. If jobVT == lapack.SVDInPlace vt is
+// of size min(m,n)×n. If jobVT == lapack.SVDOverwrite or lapack.SVDNone, vt is
+// not used.
+//
+// work is a slice for storing temporary memory, and lwork is the usable size of
+// the slice. lwork must be at least max(5*min(m,n), 3*min(m,n)+max(m,n)).
+// If lwork == -1, instead of performing Dgesvd, the optimal work length will be
+// stored into work[0]. Dgesvd will panic if the working memory has insufficient
+// storage.
+//
+// Dgesvd returns whether the decomposition successfully completed.
+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) {
+ minmn := min(m, n)
+ checkMatrix(m, n, a, lda)
+ if jobU == lapack.SVDAll {
+ checkMatrix(m, m, u, ldu)
+ } else if jobU == lapack.SVDInPlace {
+ checkMatrix(m, minmn, u, ldu)
+ }
+ if jobVT == lapack.SVDAll {
+ checkMatrix(n, n, vt, ldvt)
+ } else if jobVT == lapack.SVDInPlace {
+ checkMatrix(minmn, n, vt, ldvt)
+ }
+ if jobU == lapack.SVDOverwrite && jobVT == lapack.SVDOverwrite {
+ panic("lapack: both jobU and jobVT are lapack.SVDOverwrite")
+ }
+ if len(s) < minmn {
+ panic(badS)
+ }
+ if jobU == lapack.SVDOverwrite || jobVT == lapack.SVDOverwrite {
+ panic(noSVDO)
+ }
+ if m == 0 || n == 0 {
+ return true
+ }
+
+ wantua := jobU == lapack.SVDAll
+ wantus := jobU == lapack.SVDInPlace
+ wantuas := wantua || wantus
+ wantuo := jobU == lapack.SVDOverwrite
+ wantun := jobU == lapack.None
+
+ wantva := jobVT == lapack.SVDAll
+ wantvs := jobVT == lapack.SVDInPlace
+ wantvas := wantva || wantvs
+ wantvo := jobVT == lapack.SVDOverwrite
+ wantvn := jobVT == lapack.None
+
+ bi := blas64.Implementation()
+ var mnthr int
+
+ // Compute optimal space for subroutines.
+ maxwrk := 1
+ opts := string(jobU) + string(jobVT)
+ var wrkbl, bdspac int
+ if m >= n {
+ mnthr = impl.Ilaenv(6, "DGESVD", opts, m, n, 0, 0)
+ bdspac = 5 * n
+ impl.Dgeqrf(m, n, a, lda, nil, work, -1)
+ lwork_dgeqrf := int(work[0])
+ impl.Dorgqr(m, n, n, a, lda, nil, work, -1)
+ lwork_dorgqr_n := int(work[0])
+ impl.Dorgqr(m, m, n, a, lda, nil, work, -1)
+ lwork_dorgqr_m := int(work[0])
+ impl.Dgebrd(n, n, a, lda, s, nil, nil, nil, work, -1)
+ lwork_dgebrd := int(work[0])
+ impl.Dorgbr(lapack.ApplyP, n, n, n, a, lda, nil, work, -1)
+ lwork_dorgbr_p := int(work[0])
+ impl.Dorgbr(lapack.ApplyQ, n, n, n, a, lda, nil, work, -1)
+ lwork_dorgbr_q := int(work[0])
+
+ if m >= mnthr {
+ // m >> n
+ if wantun {
+ // Path 1
+ maxwrk = n + lwork_dgeqrf
+ maxwrk = max(maxwrk, 3*n+lwork_dgebrd)
+ if wantvo || wantvas {
+ maxwrk = max(maxwrk, 3*n+lwork_dorgbr_p)
+ }
+ maxwrk = max(maxwrk, bdspac)
+ } else if wantuo && wantvn {
+ // Path 2
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_n)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = max(n*n+wrkbl, n*n+m*n+n)
+ } else if wantuo && wantvs {
+ // Path 3
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_n)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = max(n*n+wrkbl, n*n+m*n+n)
+ } else if wantus && wantvn {
+ // Path 4
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_n)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = n*n + wrkbl
+ } else if wantus && wantvo {
+ // Path 5
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_n)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = 2*n*n + wrkbl
+ } else if wantus && wantvas {
+ // Path 6
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_n)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = n*n + wrkbl
+ } else if wantua && wantvn {
+ // Path 7
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_m)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = n*n + wrkbl
+ } else if wantua && wantvo {
+ // Path 8
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_m)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = 2*n*n + wrkbl
+ } else if wantua && wantvas {
+ // Path 9
+ wrkbl = n + lwork_dgeqrf
+ wrkbl = max(wrkbl, n+lwork_dorgqr_m)
+ wrkbl = max(wrkbl, 3*n+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, 3*n+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = n*n + wrkbl
+ }
+ } else {
+ // Path 10: m > n
+ impl.Dgebrd(m, n, a, lda, s, nil, nil, nil, work, -1)
+ lwork_dgebrd := int(work[0])
+ maxwrk = 3*n + lwork_dgebrd
+ if wantus || wantuo {
+ impl.Dorgbr(lapack.ApplyQ, m, n, n, a, lda, nil, work, -1)
+ lwork_dorgbr_q = int(work[0])
+ maxwrk = max(maxwrk, 3*n+lwork_dorgbr_q)
+ }
+ if wantua {
+ impl.Dorgbr(lapack.ApplyQ, m, m, n, a, lda, nil, work, -1)
+ lwork_dorgbr_q := int(work[0])
+ maxwrk = max(maxwrk, 3*n+lwork_dorgbr_q)
+ }
+ if !wantvn {
+ maxwrk = max(maxwrk, 3*n+lwork_dorgbr_p)
+ }
+ maxwrk = max(maxwrk, bdspac)
+ }
+ } else {
+ mnthr = impl.Ilaenv(6, "DGESVD", opts, m, n, 0, 0)
+
+ bdspac = 5 * m
+ impl.Dgelqf(m, n, a, lda, nil, work, -1)
+ lwork_dgelqf := int(work[0])
+ impl.Dorglq(n, n, m, nil, n, nil, work, -1)
+ lwork_dorglq_n := int(work[0])
+ impl.Dorglq(m, n, m, a, lda, nil, work, -1)
+ lwork_dorglq_m := int(work[0])
+ impl.Dgebrd(m, m, a, lda, s, nil, nil, nil, work, -1)
+ lwork_dgebrd := int(work[0])
+ impl.Dorgbr(lapack.ApplyP, m, m, m, a, n, nil, work, -1)
+ lwork_dorgbr_p := int(work[0])
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, a, n, nil, work, -1)
+ lwork_dorgbr_q := int(work[0])
+ if n >= mnthr {
+ // n >> m
+ if wantvn {
+ // Path 1t
+ maxwrk = m + lwork_dgelqf
+ maxwrk = max(maxwrk, 3*m+lwork_dgebrd)
+ if wantuo || wantuas {
+ maxwrk = max(maxwrk, 3*m+lwork_dorgbr_q)
+ }
+ maxwrk = max(maxwrk, bdspac)
+ } else if wantvo && wantun {
+ // Path 2t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_m)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = max(m*m+wrkbl, m*m+m*n+m)
+ } else if wantvo && wantuas {
+ // Path 3t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_m)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = max(m*m+wrkbl, m*m+m*n+m)
+ } else if wantvs && wantun {
+ // Path 4t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_m)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = m*m + wrkbl
+ } else if wantvs && wantuo {
+ // Path 5t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_m)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = 2*m*m + wrkbl
+ } else if wantvs && wantuas {
+ // Path 6t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_m)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = m*m + wrkbl
+ } else if wantva && wantun {
+ // Path 7t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_n)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = m*m + wrkbl
+ } else if wantva && wantuo {
+ // Path 8t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_n)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = 2*m*m + wrkbl
+ } else if wantva && wantuas {
+ // Path 9t
+ wrkbl = m + lwork_dgelqf
+ wrkbl = max(wrkbl, m+lwork_dorglq_n)
+ wrkbl = max(wrkbl, 3*m+lwork_dgebrd)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_p)
+ wrkbl = max(wrkbl, 3*m+lwork_dorgbr_q)
+ wrkbl = max(wrkbl, bdspac)
+ maxwrk = m*m + wrkbl
+ }
+ } else {
+ // Path 10t, n > m
+ impl.Dgebrd(m, n, a, lda, s, nil, nil, nil, work, -1)
+ lwork_dgebrd = int(work[0])
+ maxwrk = 3*m + lwork_dgebrd
+ if wantvs || wantvo {
+ impl.Dorgbr(lapack.ApplyP, m, n, m, a, n, nil, work, -1)
+ lwork_dorgbr_p = int(work[0])
+ maxwrk = max(maxwrk, 3*m+lwork_dorgbr_p)
+ }
+ if wantva {
+ impl.Dorgbr(lapack.ApplyP, n, n, m, a, n, nil, work, -1)
+ lwork_dorgbr_p = int(work[0])
+ maxwrk = max(maxwrk, 3*m+lwork_dorgbr_p)
+ }
+ if !wantun {
+ maxwrk = max(maxwrk, 3*m+lwork_dorgbr_q)
+ }
+ maxwrk = max(maxwrk, bdspac)
+ }
+ }
+
+ minWork := max(1, 5*minmn)
+ if !((wantun && m >= mnthr) || (wantvn && n >= mnthr)) {
+ minWork = max(minWork, 3*minmn+max(m, n))
+ }
+
+ if lwork != -1 {
+ if len(work) < lwork {
+ panic(badWork)
+ }
+ if lwork < minWork {
+ panic(badWork)
+ }
+ }
+ if m == 0 || n == 0 {
+ return true
+ }
+
+ maxwrk = max(maxwrk, minWork)
+ work[0] = float64(maxwrk)
+ if lwork == -1 {
+ return true
+ }
+
+ // Perform decomposition.
+ eps := dlamchE
+ smlnum := math.Sqrt(dlamchS) / eps
+ bignum := 1 / smlnum
+
+ // Scale A if max element outside range [smlnum, bignum].
+ anrm := impl.Dlange(lapack.MaxAbs, m, n, a, lda, nil)
+ var iscl bool
+ if anrm > 0 && anrm < smlnum {
+ iscl = true
+ impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, m, n, a, lda)
+ } else if anrm > bignum {
+ iscl = true
+ impl.Dlascl(lapack.General, 0, 0, anrm, bignum, m, n, a, lda)
+ }
+
+ var ie int
+ if m >= n {
+ // If A has sufficiently more rows than columns, use the QR decomposition.
+ if m >= mnthr {
+ // m >> n
+ if wantun {
+ // Path 1.
+ itau := 0
+ iwork := itau + n
+
+ // Compute A = Q * R.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Zero out below R.
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
+ ie = 0
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+ // Bidiagonalize R in A.
+ impl.Dgebrd(n, n, a, lda, s, work[ie:], work[itauq:],
+ work[itaup:], work[iwork:], lwork-iwork)
+ ncvt := 0
+ if wantvo || wantvas {
+ // Generate P^T.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, a, lda, work[itaup:],
+ work[iwork:], lwork-iwork)
+ ncvt = n
+ }
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration computing right singular vectors
+ // of A in A if desired.
+ ok = impl.Dbdsqr(blas.Upper, n, ncvt, 0, 0, s, work[ie:],
+ a, lda, work, 1, work, 1, work[iwork:])
+
+ // If right singular vectors desired in VT, copy them there.
+ if wantvas {
+ impl.Dlacpy(blas.All, n, n, a, lda, vt, ldvt)
+ }
+ } else if wantuo && wantvn {
+ // Path 2
+ panic(noSVDO)
+ } else if wantuo && wantvas {
+ // Path 3
+ panic(noSVDO)
+ } else if wantus {
+ if wantvn {
+ // Path 4
+ if lwork >= n*n+max(4*n, bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ ir := 0
+ var ldworkr int
+ if lwork >= wrkbl+lda*n {
+ ldworkr = lda
+ } else {
+ ldworkr = n
+ }
+ itau := ir + ldworkr*n
+ iwork := itau + n
+ // Compute A = Q * R.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy R to work[ir:], zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, work[ir:], ldworkr)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[ir+ldworkr:], ldworkr)
+
+ // Generate Q in A.
+ impl.Dorgqr(m, n, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in work[ir:].
+ impl.Dgebrd(n, n, work[ir:], ldworkr, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate left vectors bidiagonalizing R in work[ir:].
+ impl.Dorgbr(lapack.ApplyQ, n, n, n, work[ir:], ldworkr,
+ work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, compuing left singular
+ // vectors of R in work[ir:].
+ ok = impl.Dbdsqr(blas.Upper, n, 0, n, 0, s, work[ie:], work, 1,
+ work[ir:], ldworkr, work, 1, work[iwork:])
+
+ // Multiply Q in A by left singular vectors of R in
+ // work[ir:], storing result in U.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, a, lda,
+ work[ir:], ldworkr, 0, u, ldu)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + n
+
+ // Compute A = Q*R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, n, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Zero out below R in A.
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
+
+ // Bidiagonalize R in A.
+ impl.Dgebrd(n, n, a, lda, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply Q in U by left vectors bidiagonalizing R.
+ impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
+ a, lda, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left
+ // singular vectors of A in U.
+ ok = impl.Dbdsqr(blas.Upper, n, 0, m, 0, s, work[ie:], work, 1,
+ u, ldu, work, 1, work[iwork:])
+ }
+ } else if wantvo {
+ // Path 5
+ panic(noSVDO)
+ } else if wantvas {
+ // Path 6
+ if lwork >= n*n+max(4*n, bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ iu := 0
+ var ldworku int
+ if lwork >= wrkbl+lda*n {
+ ldworku = lda
+ } else {
+ ldworku = n
+ }
+ itau := iu + ldworku*n
+ iwork := itau + n
+
+ // Compute A = Q * R.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ // Copy R to work[iu:], zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, work[iu:], ldworku)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[iu+ldworku:], ldworku)
+
+ // Generate Q in A.
+ impl.Dorgqr(m, n, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in work[iu:], copying result to VT.
+ impl.Dgebrd(n, n, work[iu:], ldworku, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, n, n, work[iu:], ldworku, vt, ldvt)
+
+ // Generate left bidiagonalizing vectors in work[iu:].
+ impl.Dorgbr(lapack.ApplyQ, n, n, n, work[iu:], ldworku,
+ work[itauq:], work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonalizing vectors in VT.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of R in work[iu:], and computing right singular
+ // vectors of R in VT.
+ ok = impl.Dbdsqr(blas.Upper, n, n, n, 0, s, work[ie:],
+ vt, ldvt, work[iu:], ldworku, work, 1, work[iwork:])
+
+ // Multiply Q in A by left singular vectors of R in
+ // work[iu:], storing result in U.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, a, lda,
+ work[iu:], ldworku, 0, u, ldu)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + n
+
+ // Compute A = Q * R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, n, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy R to VT, zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, vt[ldvt:], ldvt)
+
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in VT.
+ impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply Q in U by left bidiagonalizing vectors in VT.
+ impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
+ vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonalizing vectors in VT.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of A in U and computing right singular vectors
+ // of A in VT.
+ ok = impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ }
+ }
+ } else if wantua {
+ if wantvn {
+ // Path 7
+ if lwork >= n*n+max(max(n+m, 4*n), bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ ir := 0
+ var ldworkr int
+ if lwork >= wrkbl+lda*n {
+ ldworkr = lda
+ } else {
+ ldworkr = n
+ }
+ itau := ir + ldworkr*n
+ iwork := itau + n
+
+ // Compute A = Q*R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Copy R to work[ir:], zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, work[ir:], ldworkr)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[ir+ldworkr:], ldworkr)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in work[ir:].
+ impl.Dgebrd(n, n, work[ir:], ldworkr, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate left bidiagonalizing vectors in work[ir:].
+ impl.Dorgbr(lapack.ApplyQ, n, n, n, work[ir:], ldworkr,
+ work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of R in work[ir:].
+ ok = impl.Dbdsqr(blas.Upper, n, 0, n, 0, s, work[ie:], work, 1,
+ work[ir:], ldworkr, work, 1, work[iwork:])
+
+ // Multiply Q in U by left singular vectors of R in
+ // work[ir:], storing result in A.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1, u, ldu,
+ work[ir:], ldworkr, 0, a, lda)
+
+ // Copy left singular vectors of A from A to U.
+ impl.Dlacpy(blas.All, m, n, a, lda, u, ldu)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + n
+
+ // Compute A = Q*R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Zero out below R in A.
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda)
+
+ // Bidiagonalize R in A.
+ impl.Dgebrd(n, n, a, lda, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply Q in U by left bidiagonalizing vectors in A.
+ impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans, m, n, n,
+ a, lda, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left
+ // singular vectors of A in U.
+ ok = impl.Dbdsqr(blas.Upper, n, 0, m, 0, s, work[ie:],
+ work, 1, u, ldu, work, 1, work[iwork:])
+ }
+ } else if wantvo {
+ // Path 8.
+ panic(noSVDO)
+ } else if wantvas {
+ // Path 9.
+ if lwork >= n*n+max(max(n+m, 4*n), bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ iu := 0
+ var ldworku int
+ if lwork >= wrkbl+lda*n {
+ ldworku = lda
+ } else {
+ ldworku = n
+ }
+ itau := iu + ldworku*n
+ iwork := itau + n
+
+ // Compute A = Q * R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy R to work[iu:], zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, work[iu:], ldworku)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, work[iu+ldworku:], ldworku)
+
+ ie = itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in work[iu:], copying result to VT.
+ impl.Dgebrd(n, n, work[iu:], ldworku, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, n, n, work[iu:], ldworku, vt, ldvt)
+
+ // Generate left bidiagonalizing vectors in work[iu:].
+ impl.Dorgbr(lapack.ApplyQ, n, n, n, work[iu:], ldworku,
+ work[itauq:], work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonalizing vectors in VT.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of R in work[iu:] and computing right
+ // singular vectors of R in VT.
+ ok = impl.Dbdsqr(blas.Upper, n, n, n, 0, s, work[ie:],
+ vt, ldvt, work[iu:], ldworku, work, 1, work[iwork:])
+
+ // Multiply Q in U by left singular vectors of R in
+ // work[iu:], storing result in A.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, n, 1,
+ u, ldu, work[iu:], ldworku, 0, a, lda)
+
+ // Copy left singular vectors of A from A to U.
+ impl.Dlacpy(blas.All, m, n, a, lda, u, ldu)
+
+ /*
+ // Bidiagonalize R in VT.
+ impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply Q in U by left bidiagonalizing vectors in VT.
+ impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans,
+ m, n, n, vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonalizing vectors in VT.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of A in U and computing right singular vectors
+ // of A in VT.
+ ok = impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ */
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + n
+
+ // Compute A = Q*R, copying result to U.
+ impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+
+ // Generate Q in U.
+ impl.Dorgqr(m, m, n, u, ldu, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy R from A to VT, zeroing out below it.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
+ impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, vt[ldvt:], ldvt)
+
+ ie := itau
+ itauq := ie + n
+ itaup := itauq + n
+ iwork = itaup + n
+
+ // Bidiagonalize R in VT.
+ impl.Dgebrd(n, n, vt, ldvt, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply Q in U by left bidiagonalizing vectors in VT.
+ impl.Dormbr(lapack.ApplyQ, blas.Right, blas.NoTrans,
+ m, n, n, vt, ldvt, work[itauq:], u, ldu, work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonizing vectors in VT.
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + n
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of A in U and computing right singular vectors
+ // of A in VT.
+ impl.Dbdsqr(blas.Upper, n, n, m, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ }
+ }
+ }
+ } else {
+ // Path 10.
+ // M at least N, but not much larger.
+ ie = 0
+ itauq := ie + n
+ itaup := itauq + n
+ iwork := itaup + n
+
+ // Bidiagonalize A.
+ impl.Dgebrd(m, n, a, lda, s, work[ie:], work[itauq:],
+ work[itaup:], work[iwork:], lwork-iwork)
+ if wantuas {
+ // Left singular vectors are desired in U. Copy result to U and
+ // generate left biadiagonalizing vectors in U.
+ impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu)
+ var ncu int
+ if wantus {
+ ncu = n
+ }
+ if wantua {
+ ncu = m
+ }
+ impl.Dorgbr(lapack.ApplyQ, m, ncu, n, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ }
+ if wantvas {
+ // Right singular vectors are desired in VT. Copy result to VT and
+ // generate left biadiagonalizing vectors in VT.
+ impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt)
+ impl.Dorgbr(lapack.ApplyP, n, n, n, vt, ldvt, work[itaup:], work[iwork:], lwork-iwork)
+ }
+ if wantuo {
+ panic(noSVDO)
+ }
+ if wantvo {
+ panic(noSVDO)
+ }
+ iwork = ie + n
+ var nru, ncvt int
+ if wantuas || wantuo {
+ nru = m
+ }
+ if wantun {
+ nru = 0
+ }
+ if wantvas || wantvo {
+ ncvt = n
+ }
+ if wantvn {
+ ncvt = 0
+ }
+ if !wantuo && !wantvo {
+ // Perform bidiagonal QR iteration, if desired, computing left
+ // singular vectors in U and right singular vectors in VT.
+ ok = impl.Dbdsqr(blas.Upper, n, ncvt, nru, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ } else {
+ // There will be two branches when the implementation is complete.
+ panic(noSVDO)
+ }
+ }
+ } else {
+ // A has more columns than rows. If A has sufficiently more columns than
+ // rows, first reduce using the LQ decomposition.
+ if n >= mnthr {
+ // n >> m.
+ if wantvn {
+ // Path 1t.
+ itau := 0
+ iwork := itau + m
+
+ // Compute A = L*Q.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Zero out above L.
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
+ ie := 0
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in A.
+ impl.Dgebrd(m, m, a, lda, s, work[ie:itauq],
+ work[itauq:itaup], work[itaup:iwork], work[iwork:], lwork-iwork)
+ if wantuo || wantuas {
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, a, lda,
+ work[itauq:], work[iwork:], lwork-iwork)
+ }
+ iwork = ie + m
+ nru := 0
+ if wantuo || wantuas {
+ nru = m
+ }
+
+ // Perform bidiagonal QR iteration, computing left singular vectors
+ // of A in A if desired.
+ ok = impl.Dbdsqr(blas.Upper, m, 0, nru, 0, s, work[ie:],
+ work, 1, a, lda, work, 1, work[iwork:])
+
+ // If left singular vectors desired in U, copy them there.
+ if wantuas {
+ impl.Dlacpy(blas.All, m, m, a, lda, u, ldu)
+ }
+ } else if wantvo && wantun {
+ // Path 2t.
+ panic(noSVDO)
+ } else if wantvo && wantuas {
+ // Path 3t.
+ panic(noSVDO)
+ } else if wantvs {
+ if wantun {
+ // Path 4t.
+ if lwork >= m*m+max(4*m, bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ ir := 0
+ var ldworkr int
+ if lwork >= wrkbl+lda*m {
+ ldworkr = lda
+ } else {
+ ldworkr = m
+ }
+ itau := ir + ldworkr*m
+ iwork := itau + m
+
+ // Compute A = L*Q.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy L to work[ir:], zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, work[ir:], ldworkr)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[ir+1:], ldworkr)
+
+ // Generate Q in A.
+ impl.Dorglq(m, n, m, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in work[ir:].
+ impl.Dgebrd(m, m, work[ir:], ldworkr, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate right vectors bidiagonalizing L in work[ir:].
+ impl.Dorgbr(lapack.ApplyP, m, m, m, work[ir:], ldworkr,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing right singular
+ // vectors of L in work[ir:].
+ ok = impl.Dbdsqr(blas.Upper, m, m, 0, 0, s, work[ie:],
+ work[ir:], ldworkr, work, 1, work, 1, work[iwork:])
+
+ // Multiply right singular vectors of L in work[ir:] by
+ // Q in A, storing result in VT.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
+ work[ir:], ldworkr, a, lda, 0, vt, ldvt)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + m
+
+ // Compute A = L*Q.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy result to VT.
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Generate Q in VT.
+ impl.Dorglq(m, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Zero out above L in A.
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
+
+ // Bidiagonalize L in A.
+ impl.Dgebrd(m, m, a, lda, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply right vectors bidiagonalizing L by Q in VT.
+ impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
+ a, lda, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing right
+ // singular vectors of A in VT.
+ ok = impl.Dbdsqr(blas.Upper, m, n, 0, 0, s, work[ie:],
+ vt, ldvt, work, 1, work, 1, work[iwork:])
+ }
+ } else if wantuo {
+ // Path 5t.
+ panic(noSVDO)
+ } else if wantuas {
+ // Path 6t.
+ if lwork >= m*m+max(4*m, bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ iu := 0
+ var ldworku int
+ if lwork >= wrkbl+lda*m {
+ ldworku = lda
+ } else {
+ ldworku = m
+ }
+ itau := iu + ldworku*m
+ iwork := itau + m
+
+ // Compute A = L*Q.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy L to work[iu:], zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, work[iu:], ldworku)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[iu+1:], ldworku)
+
+ // Generate Q in A.
+ impl.Dorglq(m, n, m, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in work[iu:], copying result to U.
+ impl.Dgebrd(m, m, work[iu:], ldworku, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, m, work[iu:], ldworku, u, ldu)
+
+ // Generate right bidiagionalizing vectors in work[iu:].
+ impl.Dorgbr(lapack.ApplyP, m, m, m, work[iu:], ldworku,
+ work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate left bidiagonalizing vectors in U.
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of L in U and computing right singular vectors of
+ // L in work[iu:].
+ ok = impl.Dbdsqr(blas.Upper, m, m, m, 0, s, work[ie:],
+ work[iu:], ldworku, u, ldu, work, 1, work[iwork:])
+
+ // Multiply right singular vectors of L in work[iu:] by
+ // Q in A, storing result in VT.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
+ work[iu:], ldworku, a, lda, 0, vt, ldvt)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + m
+
+ // Compute A = L*Q, copying result to VT.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Generate Q in VT.
+ impl.Dorglq(m, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy L to U, zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, u[1:], ldu)
+
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in U.
+ impl.Dgebrd(m, m, u, ldu, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply right bidiagonalizing vectors in U by Q in VT.
+ impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
+ u, ldu, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
+
+ // Generate left bidiagonalizing vectors in U.
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of A in U and computing right singular vectors
+ // of A in VT.
+ impl.Dbdsqr(blas.Upper, m, n, m, 0, s, work[ie:], vt, ldvt,
+ u, ldu, work, 1, work[iwork:])
+ }
+ }
+ } else if wantva {
+ if wantun {
+ // Path 7t.
+ if lwork >= m*m+max(max(n+m, 4*m), bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ ir := 0
+ var ldworkr int
+ if lwork >= wrkbl+lda*m {
+ ldworkr = lda
+ } else {
+ ldworkr = m
+ }
+ itau := ir + ldworkr*m
+ iwork := itau + m
+
+ // Compute A = L*Q, copying result to VT.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Copy L to work[ir:], zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, work[ir:], ldworkr)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[ir+1:], ldworkr)
+
+ // Generate Q in VT.
+ impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in work[ir:].
+ impl.Dgebrd(m, m, work[ir:], ldworkr, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate right bidiagonalizing vectors in work[ir:].
+ impl.Dorgbr(lapack.ApplyP, m, m, m, work[ir:], ldworkr,
+ work[itaup:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing right
+ // singular vectors of L in work[ir:].
+ ok = impl.Dbdsqr(blas.Upper, m, m, 0, 0, s, work[ie:],
+ work[ir:], ldworkr, work, 1, work, 1, work[iwork:])
+
+ // Multiply right singular vectors of L in work[ir:] by
+ // Q in VT, storing result in A.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
+ work[ir:], ldworkr, vt, ldvt, 0, a, lda)
+
+ // Copy right singular vectors of A from A to VT.
+ impl.Dlacpy(blas.All, m, n, a, lda, vt, ldvt)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + m
+ // Compute A = L * Q, copying result to VT.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Generate Q in VT.
+ impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+
+ ie := itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Zero out above L in A.
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, a[1:], lda)
+
+ // Bidiagonalize L in A.
+ impl.Dgebrd(m, m, a, lda, s, work[ie:], work[itauq:],
+ work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply right bidiagonalizing vectors in A by Q in VT.
+ impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
+ a, lda, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing right singular
+ // vectors of A in VT.
+ ok = impl.Dbdsqr(blas.Upper, m, n, 0, 0, s, work[ie:],
+ vt, ldvt, work, 1, work, 1, work[iwork:])
+ }
+ } else if wantuo {
+ panic(noSVDO)
+ } else if wantuas {
+ // Path 9t.
+ if lwork >= m*m+max(max(m+n, 4*m), bdspac) {
+ // Sufficient workspace for a fast algorithm.
+ iu := 0
+
+ var ldworku int
+ if lwork >= wrkbl+lda*m {
+ ldworku = lda
+ } else {
+ ldworku = m
+ }
+ itau := iu + ldworku*m
+ iwork := itau + m
+
+ // Generate A = L * Q copying result to VT.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Generate Q in VT.
+ impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy L to work[iu:], zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, work[iu:], ldworku)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, work[iu+1:], ldworku)
+ ie = itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in work[iu:], copying result to U.
+ impl.Dgebrd(m, m, work[iu:], ldworku, s, work[ie:],
+ work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Lower, m, m, work[iu:], ldworku, u, ldu)
+
+ // Generate right bidiagonalizing vectors in work[iu:].
+ impl.Dorgbr(lapack.ApplyP, m, m, m, work[iu:], ldworku,
+ work[itaup:], work[iwork:], lwork-iwork)
+
+ // Generate left bidiagonalizing vectors in U.
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of L in U and computing right singular vectors
+ // of L in work[iu:].
+ ok = impl.Dbdsqr(blas.Upper, m, m, m, 0, s, work[ie:],
+ work[iu:], ldworku, u, ldu, work, 1, work[iwork:])
+
+ // Multiply right singular vectors of L in work[iu:]
+ // Q in VT, storing result in A.
+ bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n, m, 1,
+ work[iu:], ldworku, vt, ldvt, 0, a, lda)
+
+ // Copy right singular vectors of A from A to VT.
+ impl.Dlacpy(blas.All, m, n, a, lda, vt, ldvt)
+ } else {
+ // Insufficient workspace for a fast algorithm.
+ itau := 0
+ iwork := itau + m
+
+ // Compute A = L * Q, copying result to VT.
+ impl.Dgelqf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork)
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+
+ // Generate Q in VT.
+ impl.Dorglq(n, n, m, vt, ldvt, work[itau:], work[iwork:], lwork-iwork)
+
+ // Copy L to U, zeroing out above it.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
+ impl.Dlaset(blas.Upper, m-1, m-1, 0, 0, u[1:], ldu)
+
+ ie = itau
+ itauq := ie + m
+ itaup := itauq + m
+ iwork = itaup + m
+
+ // Bidiagonalize L in U.
+ impl.Dgebrd(m, m, u, ldu, s, work[ie:], work[itauq:],
+ work[itaup:], work[iwork:], lwork-iwork)
+
+ // Multiply right bidiagonalizing vectors in U by Q in VT.
+ impl.Dormbr(lapack.ApplyP, blas.Left, blas.Trans, m, n, m,
+ u, ldu, work[itaup:], vt, ldvt, work[iwork:], lwork-iwork)
+
+ // Generate left bidiagonalizing vectors in U.
+ impl.Dorgbr(lapack.ApplyQ, m, m, m, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ iwork = ie + m
+
+ // Perform bidiagonal QR iteration, computing left singular
+ // vectors of A in U and computing right singular vectors
+ // of A in VT.
+ ok = impl.Dbdsqr(blas.Upper, m, n, m, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ }
+ }
+ }
+ } else {
+ // Path 10t.
+ // N at least M, but not much larger.
+ ie = 0
+ itauq := ie + m
+ itaup := itauq + m
+ iwork := itaup + m
+
+ // Bidiagonalize A.
+ impl.Dgebrd(m, n, a, lda, s, work[ie:], work[itauq:], work[itaup:], work[iwork:], lwork-iwork)
+ if wantuas {
+ // If left singular vectors desired in U, copy result to U and
+ // generate left bidiagonalizing vectors in U.
+ impl.Dlacpy(blas.Lower, m, m, a, lda, u, ldu)
+ impl.Dorgbr(lapack.ApplyQ, m, m, n, u, ldu, work[itauq:], work[iwork:], lwork-iwork)
+ }
+ if wantvas {
+ // If right singular vectors desired in VT, copy result to VT
+ // and generate right bidiagonalizing vectors in VT.
+ impl.Dlacpy(blas.Upper, m, n, a, lda, vt, ldvt)
+ var nrvt int
+ if wantva {
+ nrvt = n
+ } else {
+ nrvt = m
+ }
+ impl.Dorgbr(lapack.ApplyP, nrvt, n, m, vt, ldvt, work[itaup:], work[iwork:], lwork-iwork)
+ }
+ if wantuo {
+ panic(noSVDO)
+ }
+ if wantvo {
+ panic(noSVDO)
+ }
+ iwork = ie + m
+ var nru, ncvt int
+ if wantuas || wantuo {
+ nru = m
+ }
+ if wantvas || wantvo {
+ ncvt = n
+ }
+ if !wantuo && !wantvo {
+ // Perform bidiagonal QR iteration, if desired, computing left
+ // singular vectors in U and computing right singular vectors in
+ // VT.
+ ok = impl.Dbdsqr(blas.Lower, m, ncvt, nru, 0, s, work[ie:],
+ vt, ldvt, u, ldu, work, 1, work[iwork:])
+ } else {
+ // There will be two branches when the implementation is complete.
+ panic(noSVDO)
+ }
+ }
+ }
+ if !ok {
+ if ie > 1 {
+ for i := 0; i < minmn-1; i++ {
+ work[i+1] = work[i+ie]
+ }
+ }
+ if ie < 1 {
+ for i := minmn - 2; i >= 0; i-- {
+ work[i+1] = work[i+ie]
+ }
+ }
+ }
+ // Undo scaling if necessary.
+ if iscl {
+ if anrm > bignum {
+ impl.Dlascl(lapack.General, 0, 0, bignum, anrm, minmn, 1, s, minmn)
+ }
+ if !ok && anrm > bignum {
+ impl.Dlascl(lapack.General, 0, 0, bignum, anrm, minmn-1, 1, work[minmn:], minmn)
+ }
+ if anrm < smlnum {
+ impl.Dlascl(lapack.General, 0, 0, smlnum, anrm, minmn, 1, s, minmn)
+ }
+ if !ok && anrm < smlnum {
+ impl.Dlascl(lapack.General, 0, 0, smlnum, anrm, minmn-1, 1, work[minmn:], minmn)
+ }
+ }
+ work[0] = float64(maxwrk)
+ return ok
+}
OSDN Git Service
