mat: change factorization inputs to use bit types (#872)

* mat: change factorization inputs to use bit types

Fixes #756 and #748.
This commit is contained in:
Brendan Tracey
2019-03-23 17:20:14 +00:00
committed by GitHub
parent 9dd6aa72d7
commit 1d8caee34e
7 changed files with 205 additions and 123 deletions

View File

@@ -13,42 +13,3 @@ const (
// Lower specifies a lower triangular matrix. // Lower specifies a lower triangular matrix.
Lower TriKind = false Lower TriKind = false
) )
// SVDKind specifies the treatment of singular vectors during an SVD
// factorization.
type SVDKind int
const (
// SVDNone specifies that no singular vectors should be computed during
// the decomposition.
SVDNone SVDKind = iota + 1
// SVDThin computes the thin singular vectors, that is, it computes
// A = U~ * Σ * V~^T
// where U~ is of size m×min(m,n), Σ is a diagonal matrix of size min(m,n)×min(m,n)
// and V~ is of size n×min(m,n).
SVDThin
// SVDFull computes the full singular value decomposition,
// A = U * Σ * V^T
// where U is of size m×m, Σ is an m×n diagonal matrix, and V is an n×n matrix.
SVDFull
)
// GSVDKind specifies the treatment of singular vectors during a GSVD
// factorization.
type GSVDKind int
const (
// GSVDU specifies that the U singular vectors should be computed during
// the decomposition.
GSVDU GSVDKind = 1 << iota
// GSVDV specifies that the V singular vectors should be computed during
// the decomposition.
GSVDV
// GSVDQ specifies that the Q singular vectors should be computed during
// the decomposition.
GSVDQ
// GSVDNone specifies that no singular vector should be computed during
// the decomposition.
GSVDNone
)

View File

@@ -104,12 +104,25 @@ func (m *Dense) EigenvectorsSym(e *EigenSym) {
m.Copy(e.vectors) m.Copy(e.vectors)
} }
// EigenKind specifies the computation of eigenvectors during factorization.
type EigenKind int
const (
// EigenNone specifies to not compute any eigenvectors.
EigenNone EigenKind = 0
// EigenLeft specifies to compute the left eigenvectors.
EigenLeft EigenKind = 1 << iota
// EigenRight specifies to compute the right eigenvectors.
EigenRight
// EigenBoth is a convenience value for computing both eigenvectors.
EigenBoth EigenKind = EigenLeft | EigenRight
)
// Eigen is a type for creating and using the eigenvalue decomposition of a dense matrix. // Eigen is a type for creating and using the eigenvalue decomposition of a dense matrix.
type Eigen struct { type Eigen struct {
n int // The size of the factorized matrix. n int // The size of the factorized matrix.
right bool // have the right eigenvectors been computed kind EigenKind
left bool // have the left eigenvectors been computed
values []complex128 values []complex128
rVectors *CDense rVectors *CDense
@@ -118,7 +131,7 @@ type Eigen struct {
// succFact returns whether the receiver contains a successful factorization. // succFact returns whether the receiver contains a successful factorization.
func (e *Eigen) succFact() bool { func (e *Eigen) succFact() bool {
return len(e.values) != 0 return e.n != 0
} }
// Factorize computes the eigenvalues of the square matrix a, and optionally // Factorize computes the eigenvalues of the square matrix a, and optionally
@@ -135,13 +148,17 @@ func (e *Eigen) succFact() bool {
// //
// Typically eigenvectors refer to right eigenvectors. // Typically eigenvectors refer to right eigenvectors.
// //
// In all cases, Factorize computes the eigenvalues of the matrix. If right and left // In all cases, Factorize computes the eigenvalues of the matrix. kind
// are true, then the right and left eigenvectors will be computed, respectively. // specifies which of the eigenvectors, if any, to compute. See the EigenKind
// documentation for more information.
// Eigen panics if the input matrix is not square. // Eigen panics if the input matrix is not square.
// //
// Factorize returns whether the decomposition succeeded. If the decomposition // Factorize returns whether the decomposition succeeded. If the decomposition
// failed, methods that require a successful factorization will panic. // failed, methods that require a successful factorization will panic.
func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) { func (e *Eigen) Factorize(a Matrix, kind EigenKind) (ok bool) {
// kill previous factorization.
e.n = 0
e.kind = 0
// Copy a because it is modified during the Lapack call. // Copy a because it is modified during the Lapack call.
r, c := a.Dims() r, c := a.Dims()
if r != c { if r != c {
@@ -150,6 +167,9 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) {
var sd Dense var sd Dense
sd.Clone(a) sd.Clone(a)
left := kind&EigenLeft != 0
right := kind&EigenRight != 0
var vl, vr Dense var vl, vr Dense
jobvl := lapack.LeftEVNone jobvl := lapack.LeftEVNone
jobvr := lapack.RightEVNone jobvr := lapack.RightEVNone
@@ -183,8 +203,7 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) {
return false return false
} }
e.n = r e.n = r
e.right = right e.kind = kind
e.left = left
// Construct complex eigenvalues from float64 data. // Construct complex eigenvalues from float64 data.
values := make([]complex128, r) values := make([]complex128, r)
@@ -212,6 +231,15 @@ func (e *Eigen) Factorize(a Matrix, left, right bool) (ok bool) {
return true return true
} }
// Kind returns the EigenKind of the decomposition. If no decomposition has been
// computed, Kind returns -1.
func (e *Eigen) Kind() EigenKind {
if !e.succFact() {
return -1
}
return e.kind
}
// Values extracts the eigenvalues of the factorized matrix. If dst is // Values extracts the eigenvalues of the factorized matrix. If dst is
// non-nil, the values are stored in-place into dst. In this case // non-nil, the values are stored in-place into dst. In this case
// dst must have length n, otherwise Values will panic. If dst is // dst must have length n, otherwise Values will panic. If dst is
@@ -281,7 +309,7 @@ func (e *Eigen) VectorsTo(dst *CDense) *CDense {
if !e.succFact() { if !e.succFact() {
panic(badFact) panic(badFact)
} }
if !e.right { if e.kind&EigenRight == 0 {
panic(badNoVect) panic(badNoVect)
} }
if dst == nil { if dst == nil {
@@ -303,7 +331,7 @@ func (e *Eigen) LeftVectorsTo(dst *CDense) *CDense {
if !e.succFact() { if !e.succFact() {
panic(badFact) panic(badFact)
} }
if !e.left { if e.kind&EigenLeft == 0 {
panic(badNoVect) panic(badNoVect)
} }
if dst == nil { if dst == nil {

View File

@@ -50,7 +50,7 @@ func ExampleEigen() {
fmt.Printf("A = %v\n\n", mat.Formatted(a, mat.Prefix(" "))) fmt.Printf("A = %v\n\n", mat.Formatted(a, mat.Prefix(" ")))
var eig mat.Eigen var eig mat.Eigen
ok := eig.Factorize(a, true, false) ok := eig.Factorize(a, mat.EigenLeft)
if !ok { if !ok {
log.Fatal("Eigendecomposition failed") log.Fatal("Eigendecomposition failed")
} }

View File

@@ -63,13 +63,13 @@ func TestEigen(t *testing.T) {
}, },
} { } {
var e1, e2, e3, e4 Eigen var e1, e2, e3, e4 Eigen
ok := e1.Factorize(test.a, true, true) ok := e1.Factorize(test.a, EigenBoth)
if !ok { if !ok {
panic("bad factorization") panic("bad factorization")
} }
e2.Factorize(test.a, false, true) e2.Factorize(test.a, EigenRight)
e3.Factorize(test.a, true, false) e3.Factorize(test.a, EigenLeft)
e4.Factorize(test.a, false, false) e4.Factorize(test.a, EigenNone)
v1 := e1.Values(nil) v1 := e1.Values(nil)
if !cmplxEqualTol(v1, test.values, 1e-14) { if !cmplxEqualTol(v1, test.values, 1e-14) {

View File

@@ -11,6 +11,29 @@ import (
"gonum.org/v1/gonum/lapack/lapack64" "gonum.org/v1/gonum/lapack/lapack64"
) )
// GSVDKind specifies the treatment of singular vectors during a GSVD
// factorization.
type GSVDKind int
const (
// GSVDNone specifies that no singular vectors should be computed during
// the decomposition.
GSVDNone GSVDKind = 0
// GSVDU specifies that the U singular vectors should be computed during
// the decomposition.
GSVDU GSVDKind = 1 << iota
// GSVDV specifies that the V singular vectors should be computed during
// the decomposition.
GSVDV
// GSVDQ specifies that the Q singular vectors should be computed during
// the decomposition.
GSVDQ
// GSVDAll is a convenience value for computing all of the singular vectors.
GSVDAll = GSVDU | GSVDV | GSVDQ
)
// GSVD is a type for creating and using the Generalized Singular Value Decomposition // GSVD is a type for creating and using the Generalized Singular Value Decomposition
// (GSVD) of a matrix. // (GSVD) of a matrix.
// //
@@ -28,12 +51,17 @@ type GSVD struct {
iwork []int iwork []int
} }
// succFact returns whether the receiver contains a successful factorization.
func (gsvd *GSVD) succFact() bool {
return gsvd.r != 0
}
// Factorize computes the generalized singular value decomposition (GSVD) of the input // Factorize computes the generalized singular value decomposition (GSVD) of the input
// the r×c matrix A and the p×c matrix B. The singular values of A and B are computed // the r×c matrix A and the p×c matrix B. The singular values of A and B are computed
// in all cases, while the singular vectors are optionally computed depending on the // in all cases, while the singular vectors are optionally computed depending on the
// input kind. // input kind.
// //
// The full singular value decomposition (kind == GSVDU|GSVDV|GSVDQ) deconstructs A and B as // The full singular value decomposition (kind == GSVDAll) deconstructs A and B as
// A = U * Σ₁ * [ 0 R ] * Q^T // A = U * Σ₁ * [ 0 R ] * Q^T
// //
// B = V * Σ₂ * [ 0 R ] * Q^T // B = V * Σ₂ * [ 0 R ] * Q^T
@@ -49,6 +77,10 @@ type GSVD struct {
// Factorize returns whether the decomposition succeeded. If the decomposition // Factorize returns whether the decomposition succeeded. If the decomposition
// failed, routines that require a successful factorization will panic. // failed, routines that require a successful factorization will panic.
func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) { func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) {
// kill the previous decomposition
gsvd.r = 0
gsvd.kind = 0
r, c := a.Dims() r, c := a.Dims()
gsvd.r, gsvd.c = r, c gsvd.r, gsvd.c = r, c
p, c := b.Dims() p, c := b.Dims()
@@ -64,7 +96,7 @@ func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) {
jobU = lapack.GSVDNone jobU = lapack.GSVDNone
jobV = lapack.GSVDNone jobV = lapack.GSVDNone
jobQ = lapack.GSVDNone jobQ = lapack.GSVDNone
case (GSVDU|GSVDV|GSVDQ)&kind != 0: case GSVDAll&kind != 0:
if GSVDU&kind != 0 { if GSVDU&kind != 0 {
jobU = lapack.GSVDU jobU = lapack.GSVDU
gsvd.u = blas64.General{ gsvd.u = blas64.General{
@@ -115,9 +147,12 @@ func (gsvd *GSVD) Factorize(a, b Matrix, kind GSVDKind) (ok bool) {
return ok return ok
} }
// Kind returns the matrix.GSVDKind of the decomposition. If no decomposition has been // Kind returns the GSVDKind of the decomposition. If no decomposition has been
// computed, Kind returns 0. // computed, Kind returns -1.
func (gsvd *GSVD) Kind() GSVDKind { func (gsvd *GSVD) Kind() GSVDKind {
if !gsvd.succFact() {
return -1
}
return gsvd.kind return gsvd.kind
} }
@@ -134,8 +169,8 @@ func (gsvd *GSVD) Rank() (k, l int) {
// //
// GeneralizedValues will panic if the receiver does not contain a successful factorization. // GeneralizedValues will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 { func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
c := gsvd.c c := gsvd.c
@@ -159,8 +194,8 @@ func (gsvd *GSVD) GeneralizedValues(v []float64) []float64 {
// //
// ValuesA will panic if the receiver does not contain a successful factorization. // ValuesA will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) ValuesA(s []float64) []float64 { func (gsvd *GSVD) ValuesA(s []float64) []float64 {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
c := gsvd.c c := gsvd.c
@@ -184,8 +219,8 @@ func (gsvd *GSVD) ValuesA(s []float64) []float64 {
// //
// ValuesB will panic if the receiver does not contain a successful factorization. // ValuesB will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) ValuesB(s []float64) []float64 { func (gsvd *GSVD) ValuesB(s []float64) []float64 {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
c := gsvd.c c := gsvd.c
@@ -207,8 +242,8 @@ func (gsvd *GSVD) ValuesB(s []float64) []float64 {
// //
// ZeroRTo will panic if the receiver does not contain a successful factorization. // ZeroRTo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense { func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
c := gsvd.c c := gsvd.c
@@ -245,8 +280,8 @@ func (gsvd *GSVD) ZeroRTo(dst *Dense) *Dense {
// //
// SigmaATo will panic if the receiver does not contain a successful factorization. // SigmaATo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense { func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
k := gsvd.k k := gsvd.k
@@ -271,8 +306,8 @@ func (gsvd *GSVD) SigmaATo(dst *Dense) *Dense {
// //
// SigmaBTo will panic if the receiver does not contain a successful factorization. // SigmaBTo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense { func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense {
if gsvd.kind == 0 { if !gsvd.succFact() {
panic("gsvd: no decomposition computed") panic(badFact)
} }
r := gsvd.r r := gsvd.r
p := gsvd.p p := gsvd.p
@@ -298,6 +333,9 @@ func (gsvd *GSVD) SigmaBTo(dst *Dense) *Dense {
// //
// UTo will panic if the receiver does not contain a successful factorization. // UTo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) UTo(dst *Dense) *Dense { func (gsvd *GSVD) UTo(dst *Dense) *Dense {
if !gsvd.succFact() {
panic(badFact)
}
if gsvd.kind&GSVDU == 0 { if gsvd.kind&GSVDU == 0 {
panic("mat: improper GSVD kind") panic("mat: improper GSVD kind")
} }
@@ -324,6 +362,9 @@ func (gsvd *GSVD) UTo(dst *Dense) *Dense {
// //
// VTo will panic if the receiver does not contain a successful factorization. // VTo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) VTo(dst *Dense) *Dense { func (gsvd *GSVD) VTo(dst *Dense) *Dense {
if !gsvd.succFact() {
panic(badFact)
}
if gsvd.kind&GSVDV == 0 { if gsvd.kind&GSVDV == 0 {
panic("mat: improper GSVD kind") panic("mat: improper GSVD kind")
} }
@@ -350,6 +391,9 @@ func (gsvd *GSVD) VTo(dst *Dense) *Dense {
// //
// QTo will panic if the receiver does not contain a successful factorization. // QTo will panic if the receiver does not contain a successful factorization.
func (gsvd *GSVD) QTo(dst *Dense) *Dense { func (gsvd *GSVD) QTo(dst *Dense) *Dense {
if !gsvd.succFact() {
panic(badFact)
}
if gsvd.kind&GSVDQ == 0 { if gsvd.kind&GSVDQ == 0 {
panic("mat: improper GSVD kind") panic("mat: improper GSVD kind")
} }

View File

@@ -24,6 +24,11 @@ type HOGSVD struct {
err error err error
} }
// succFact returns whether the receiver contains a successful factorization.
func (gsvd *HOGSVD) succFact() bool {
return gsvd.n != 0
}
// Factorize computes the higher order generalized singular value decomposition (HOGSVD) // Factorize computes the higher order generalized singular value decomposition (HOGSVD)
// of the n input r_i×c column tall matrices in m. HOGSV extends the GSVD case from 2 to n // of the n input r_i×c column tall matrices in m. HOGSV extends the GSVD case from 2 to n
// input matrices. // input matrices.
@@ -96,7 +101,7 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
s.Scale(1/float64(len(m)*(len(m)-1)), s) s.Scale(1/float64(len(m)*(len(m)-1)), s)
var eig Eigen var eig Eigen
ok = eig.Factorize(s.T(), false, true) ok = eig.Factorize(s.T(), EigenRight)
if !ok { if !ok {
gsvd.err = errors.New("hogsvd: eigen decomposition failed") gsvd.err = errors.New("hogsvd: eigen decomposition failed")
return false return false
@@ -157,8 +162,8 @@ func (gsvd *HOGSVD) Len() int {
// //
// UTo will panic if the receiver does not contain a successful factorization. // UTo will panic if the receiver does not contain a successful factorization.
func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense { func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
if gsvd.n == 0 { if !gsvd.succFact() {
panic("hogsvd: unsuccessful factorization") panic(badFact)
} }
if n < 0 || gsvd.n <= n { if n < 0 || gsvd.n <= n {
panic("hogsvd: invalid index") panic("hogsvd: invalid index")
@@ -187,8 +192,8 @@ func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense {
// //
// Values will panic if the receiver does not contain a successful factorization. // Values will panic if the receiver does not contain a successful factorization.
func (gsvd *HOGSVD) Values(s []float64, n int) []float64 { func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
if gsvd.n == 0 { if !gsvd.succFact() {
panic("hogsvd: unsuccessful factorization") panic(badFact)
} }
if n < 0 || gsvd.n <= n { if n < 0 || gsvd.n <= n {
panic("hogsvd: invalid index") panic("hogsvd: invalid index")
@@ -214,8 +219,8 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 {
// //
// VTo will panic if the receiver does not contain a successful factorization. // VTo will panic if the receiver does not contain a successful factorization.
func (gsvd *HOGSVD) VTo(dst *Dense) *Dense { func (gsvd *HOGSVD) VTo(dst *Dense) *Dense {
if gsvd.n == 0 { if !gsvd.succFact() {
panic("hogsvd: unsuccessful factorization") panic(badFact)
} }
if dst == nil { if dst == nil {
r, c := gsvd.v.Dims() r, c := gsvd.v.Dims()

View File

@@ -20,6 +20,35 @@ type SVD struct {
vt blas64.General vt blas64.General
} }
// SVDKind specifies the treatment of singular vectors during an SVD
// factorization.
type SVDKind int
const (
// SVDNone specifies that no singular vectors should be computed during
// the decomposition.
SVDNone SVDKind = 0
// SVDThinU specifies the thin decomposition for U should be computed.
SVDThinU SVDKind = 1 << (iota - 1)
// SVDFullU specifies the full decomposition for U should be computed.
SVDFullU
// SVDThinV specifies the thin decomposition for V should be computed.
SVDThinV
// SVDFullV specifies the full decomposition for V should be computed.
SVDFullV
// SVDThin is a convenience value for computing both thin vectors.
SVDThin SVDKind = SVDThinU | SVDThinV
// SVDThin is a convenience value for computing both full vectors.
SVDFull SVDKind = SVDFullU | SVDFullV
)
// succFact returns whether the receiver contains a successful factorization.
func (svd *SVD) succFact() bool {
return len(svd.s) != 0
}
// Factorize computes the singular value decomposition (SVD) of the input matrix A. // Factorize computes the singular value decomposition (SVD) of the input matrix A.
// The singular values of A are computed in all cases, while the singular // The singular values of A are computed in all cases, while the singular
// vectors are optionally computed depending on the input kind. // vectors are optionally computed depending on the input kind.
@@ -32,61 +61,67 @@ type SVD struct {
// The first min(m,n) columns of U and V are, respectively, the left and right // The first min(m,n) columns of U and V are, respectively, the left and right
// singular vectors of A. // singular vectors of A.
// //
// It is frequently not necessary to compute the full SVD. Computation time and // Significant storage space can be saved by using the thin representation of
// storage costs can be reduced using the appropriate kind. Only the singular // the SVD (kind == SVDThin) instead of the full SVD, especially if
// values can be computed (kind == SVDNone), or a "thin" representation of the // m >> n or m << n. The thin SVD finds
// orthogonal matrices U and V (kind = SVDThin). The thin representation can // A = U~ * Σ * V~^T
// save a significant amount of memory if m >> n or m << n. // where U~ is of size m×min(m,n), Σ is a diagonal matrix of size min(m,n)×min(m,n)
// and V~ is of size n×min(m,n).
// //
// Factorize returns whether the decomposition succeeded. If the decomposition // Factorize returns whether the decomposition succeeded. If the decomposition
// failed, routines that require a successful factorization will panic. // failed, routines that require a successful factorization will panic.
func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) { func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) {
// kill previous factorization
svd.s = svd.s[:0]
svd.kind = kind
m, n := a.Dims() m, n := a.Dims()
var jobU, jobVT lapack.SVDJob var jobU, jobVT lapack.SVDJob
switch kind {
default: // TODO(btracey): This code should be modified to have the smaller
panic("svd: bad input kind") // matrix written in-place into aCopy when the lapack/native/dgesvd
case SVDNone: // implementation is complete.
svd.u.Stride = 1 switch {
svd.vt.Stride = 1 case kind&SVDFullU != 0:
jobU = lapack.SVDNone jobU = lapack.SVDAll
jobVT = lapack.SVDNone
case SVDFull:
// TODO(btracey): This code should be modified to have the smaller
// matrix written in-place into aCopy when the lapack/native/dgesvd
// implementation is complete.
svd.u = blas64.General{ svd.u = blas64.General{
Rows: m, Rows: m,
Cols: m, Cols: m,
Stride: m, Stride: m,
Data: use(svd.u.Data, m*m), Data: use(svd.u.Data, m*m),
} }
svd.vt = blas64.General{ case kind&SVDThinU != 0:
Rows: n, jobU = lapack.SVDStore
Cols: n,
Stride: n,
Data: use(svd.vt.Data, n*n),
}
jobU = lapack.SVDAll
jobVT = lapack.SVDAll
case SVDThin:
// TODO(btracey): This code should be modified to have the larger
// matrix written in-place into aCopy when the lapack/native/dgesvd
// implementation is complete.
svd.u = blas64.General{ svd.u = blas64.General{
Rows: m, Rows: m,
Cols: min(m, n), Cols: min(m, n),
Stride: min(m, n), Stride: min(m, n),
Data: use(svd.u.Data, m*min(m, n)), Data: use(svd.u.Data, m*min(m, n)),
} }
default:
svd.u.Stride = 1
jobU = lapack.SVDNone
}
switch {
case kind&SVDFullV != 0:
svd.vt = blas64.General{
Rows: n,
Cols: n,
Stride: n,
Data: use(svd.vt.Data, n*n),
}
jobVT = lapack.SVDAll
case kind&SVDThinV != 0:
svd.vt = blas64.General{ svd.vt = blas64.General{
Rows: min(m, n), Rows: min(m, n),
Cols: n, Cols: n,
Stride: n, Stride: n,
Data: use(svd.vt.Data, min(m, n)*n), Data: use(svd.vt.Data, min(m, n)*n),
} }
jobU = lapack.SVDStore
jobVT = lapack.SVDStore jobVT = lapack.SVDStore
default:
svd.vt.Stride = 1
jobVT = lapack.SVDNone
} }
// A is destroyed on call, so copy the matrix. // A is destroyed on call, so copy the matrix.
@@ -105,17 +140,20 @@ func (svd *SVD) Factorize(a Matrix, kind SVDKind) (ok bool) {
return ok return ok
} }
// Kind returns the matrix.SVDKind of the decomposition. If no decomposition has been // Kind returns the SVDKind of the decomposition. If no decomposition has been
// computed, Kind returns 0. // computed, Kind returns -1.
func (svd *SVD) Kind() SVDKind { func (svd *SVD) Kind() SVDKind {
if !svd.succFact() {
return -1
}
return svd.kind return svd.kind
} }
// Cond returns the 2-norm condition number for the factorized matrix. Cond will // Cond returns the 2-norm condition number for the factorized matrix. Cond will
// panic if the receiver does not contain a successful factorization. // panic if the receiver does not contain a successful factorization.
func (svd *SVD) Cond() float64 { func (svd *SVD) Cond() float64 {
if svd.kind == 0 { if !svd.succFact() {
panic("svd: no decomposition computed") panic(badFact)
} }
return svd.s[0] / svd.s[len(svd.s)-1] return svd.s[0] / svd.s[len(svd.s)-1]
} }
@@ -129,8 +167,8 @@ func (svd *SVD) Cond() float64 {
// //
// Values will panic if the receiver does not contain a successful factorization. // Values will panic if the receiver does not contain a successful factorization.
func (svd *SVD) Values(s []float64) []float64 { func (svd *SVD) Values(s []float64) []float64 {
if svd.kind == 0 { if !svd.succFact() {
panic("svd: no decomposition computed") panic(badFact)
} }
if s == nil { if s == nil {
s = make([]float64, len(svd.s)) s = make([]float64, len(svd.s))
@@ -147,13 +185,16 @@ func (svd *SVD) Values(s []float64) []float64 {
// values as returned from SVD.Values. // values as returned from SVD.Values.
// //
// If dst is not nil, U is stored in-place into dst, and dst must have size // If dst is not nil, U is stored in-place into dst, and dst must have size
// m×m if svd.Kind() == SVDFull, size m×min(m,n) if svd.Kind() == SVDThin, and // m×m if the full U was computed, size m×min(m,n) if the thin U was computed,
// UTo panics otherwise. If dst is nil, a new matrix of the appropriate size is // and UTo panics otherwise. If dst is nil, a new matrix of the appropriate size
// allocated and returned. // is allocated and returned.
func (svd *SVD) UTo(dst *Dense) *Dense { func (svd *SVD) UTo(dst *Dense) *Dense {
if !svd.succFact() {
panic(badFact)
}
kind := svd.kind kind := svd.kind
if kind != SVDFull && kind != SVDThin { if kind&SVDThinU == 0 && kind&SVDFullU == 0 {
panic("mat: improper SVD kind") panic("svd: u not computed during factorization")
} }
r := svd.u.Rows r := svd.u.Rows
c := svd.u.Cols c := svd.u.Cols
@@ -178,13 +219,16 @@ func (svd *SVD) UTo(dst *Dense) *Dense {
// values as returned from SVD.Values. // values as returned from SVD.Values.
// //
// If dst is not nil, V is stored in-place into dst, and dst must have size // If dst is not nil, V is stored in-place into dst, and dst must have size
// n×n if svd.Kind() == SVDFull, size n×min(m,n) if svd.Kind() == SVDThin, and // n×n if the full V was computed, size n×min(m,n) if the thin V was computed,
// VTo panics otherwise. If dst is nil, a new matrix of the appropriate size is // and VTo panics otherwise. If dst is nil, a new matrix of the appropriate size
// allocated and returned. // is allocated and returned.
func (svd *SVD) VTo(dst *Dense) *Dense { func (svd *SVD) VTo(dst *Dense) *Dense {
if !svd.succFact() {
panic(badFact)
}
kind := svd.kind kind := svd.kind
if kind != SVDFull && kind != SVDThin { if kind&SVDThinU == 0 && kind&SVDFullV == 0 {
panic("mat: improper SVD kind") panic("svd: v not computed during factorization")
} }
r := svd.vt.Rows r := svd.vt.Rows
c := svd.vt.Cols c := svd.vt.Cols