Adddgesvd

cgo/lapack.go

@@ -24,6 +24,7 @@ const (
 	badLdA          = "lapack: index of a out of range"
 	badNorm         = "lapack: bad norm"
 	badPivot        = "lapack: bad pivot"
+	badS            = "lapack: s has insufficient length"
 	badSide         = "lapack: bad side"
 	badSlice        = "lapack: bad input slice length"
 	badStore        = "lapack: bad store"
@@ -497,6 +498,85 @@ func (impl Implementation) Dgels(trans blas.Transpose, m, n, nrhs int, a []float
 	return clapack.Dgels(trans, m, n, nrhs, a, lda, b, ldb)
 }
 
+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 columnwise. 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 rowwise. 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.
+//
+// The C interface does not support providing temporary storage. To provide compatibility
+// with native, lwork == -1 will not run Dgesvd but will instead write the minimum
+// work necessary to work[0]. If len(work) < lwork, Dgeqrf will panic.
+//
+// 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) {
+	checkMatrix(m, n, a, lda)
+	if jobU == lapack.SVDAll {
+		checkMatrix(m, m, u, ldu)
+	} else if jobU == lapack.SVDInPlace {
+		checkMatrix(m, min(m, n), u, ldu)
+	}
+	if jobVT == lapack.SVDAll {
+		checkMatrix(n, n, vt, ldvt)
+	} else if jobVT == lapack.SVDInPlace {
+		checkMatrix(min(m, n), n, vt, ldvt)
+	}
+	if jobU == lapack.SVDOverwrite && jobVT == lapack.SVDOverwrite {
+		panic(noSVDO)
+	}
+	if len(s) < min(m, n) {
+		panic(badS)
+	}
+	if jobU == lapack.SVDOverwrite || jobVT == lapack.SVDOverwrite {
+		panic("lapack: SVD not coded to overwrite original matrix")
+	}
+	minWork := max(5*min(m, n), 3*min(m, n)+max(m, n))
+	if lwork != -1 {
+		if len(work) < lwork {
+			panic(badWork)
+		}
+		if lwork < minWork {
+			panic(badWork)
+		}
+	}
+	if lwork == -1 {
+		work[0] = float64(minWork)
+		return true
+	}
+	return clapack.Dgesvd(lapack.Job(jobU), lapack.Job(jobVT), m, n, a, lda, s, u, ldu, vt, ldvt, work[1:])
+}
+
 // Dgetf2 computes the LU decomposition of the m×n matrix A.
 // The LU decomposition is a factorization of a into
 //  A = P * L * U

cgo/lapack_test.go

@@ -88,6 +88,10 @@ func TestDgeqrf(t *testing.T) {
 	testlapack.DgeqrfTest(t, impl)
 }
 
+func TestDgesvd(t *testing.T) {
+	testlapack.DgesvdTest(t, impl)
+}
+
 func TestDgetf2(t *testing.T) {
 	testlapack.Dgetf2Test(t, impl)
 }

lapack.go

@@ -98,3 +98,13 @@ const (
 	ApplyP DecompUpdate = 'P'
 	ApplyQ DecompUpdate = 'Q'
 )
+
+// SVDJob specifies the singular vector computation type for SVD.
+type SVDJob byte
+
+const (
+	SVDAll       SVDJob = 'A' // Compute all singular vectors
+	SVDInPlace          = 'S' // Compute the first singular vectors and store them in provided storage.
+	SVDOverwrite        = 'O' // Compute the singular vectors and store them in input matrix
+	SVDNone             = 'N' // Do not compute singular vectors
+)

native/dbdsqr.go

@@ -123,6 +123,7 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v
 		for i := 0; i < n-1; i++ {
 			smax = math.Max(smax, math.Abs(e[i]))
 		}
+
 		var sminl float64
 		var thresh float64
 		if tol >= 0 {

native/dgebrd.go

@@ -86,7 +86,7 @@ func (impl Implementation) Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, ta
 				if lwork >= (m+n)*nbmin {
 					nb = lwork / (m + n)
 				} else {
-					nb = 1
+					nb = minmn
 					nx = minmn
 				}
 			}