From 4fc257f2c0f7803846a33f5e7ebb705dfd07b445 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:35:20 +0000 Subject: [PATCH] build(deps): bump gonum.org/v1/gonum from 0.13.0 to 0.15.1 Bumps [gonum.org/v1/gonum](https://github.com/gonum/gonum) from 0.13.0 to 0.15.1. - [Release notes](https://github.com/gonum/gonum/releases) - [Commits](https://github.com/gonum/gonum/compare/v0.13.0...v0.15.1) --- updated-dependencies: - dependency-name: gonum.org/v1/gonum dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- go.mod | 4 +- go.sum | 10 +- .../github.com/google/go-cmp/cmp/compare.go | 38 +- .../cmp/{export_unsafe.go => export.go} | 5 - .../google/go-cmp/cmp/export_panic.go | 16 - .../value/{pointer_unsafe.go => pointer.go} | 3 - .../cmp/internal/value/pointer_purego.go | 34 - .../github.com/google/go-cmp/cmp/options.go | 84 +-- vendor/github.com/google/go-cmp/cmp/path.go | 46 +- .../google/go-cmp/cmp/report_reflect.go | 2 +- vendor/gonum.org/v1/gonum/AUTHORS | 12 + vendor/gonum.org/v1/gonum/CONTRIBUTORS | 12 + .../gonum.org/v1/gonum/blas/blas64/blas64.go | 2 +- vendor/gonum.org/v1/gonum/blas/blas64/conv.go | 14 - .../gonum.org/v1/gonum/blas/cblas128/conv.go | 14 - vendor/gonum.org/v1/gonum/blas/gonum/gonum.go | 14 - vendor/gonum.org/v1/gonum/floats/floats.go | 7 +- .../gonum.org/v1/gonum/lapack/gonum/dbdsqr.go | 14 +- .../gonum.org/v1/gonum/lapack/gonum/dgecon.go | 30 +- .../gonum.org/v1/gonum/lapack/gonum/dgels.go | 2 +- .../gonum.org/v1/gonum/lapack/gonum/dgeqp3.go | 43 +- .../gonum.org/v1/gonum/lapack/gonum/dgeqr2.go | 6 +- .../gonum.org/v1/gonum/lapack/gonum/dgeqrf.go | 8 +- .../gonum.org/v1/gonum/lapack/gonum/dgesvd.go | 34 +- .../gonum.org/v1/gonum/lapack/gonum/dgetf2.go | 28 +- .../gonum.org/v1/gonum/lapack/gonum/dgetrf.go | 26 +- .../gonum.org/v1/gonum/lapack/gonum/dgghrd.go | 125 ++++ .../v1/gonum/lapack/gonum/dggsvp3.go | 13 +- .../gonum.org/v1/gonum/lapack/gonum/dhseqr.go | 2 +- .../gonum.org/v1/gonum/lapack/gonum/dlahqr.go | 26 +- .../gonum.org/v1/gonum/lapack/gonum/dlanhs.go | 78 ++ .../v1/gonum/lapack/gonum/dlaqr04.go | 20 +- .../gonum.org/v1/gonum/lapack/gonum/dlaqr5.go | 700 ++++++++---------- .../gonum.org/v1/gonum/lapack/gonum/dlassq.go | 26 +- .../gonum.org/v1/gonum/lapack/gonum/dlatbs.go | 2 +- .../gonum.org/v1/gonum/lapack/gonum/dlatrs.go | 59 +- .../gonum.org/v1/gonum/lapack/gonum/dorg2r.go | 6 +- .../gonum.org/v1/gonum/lapack/gonum/dorgbr.go | 4 +- .../gonum.org/v1/gonum/lapack/gonum/dorgl2.go | 18 +- .../gonum.org/v1/gonum/lapack/gonum/dorglq.go | 26 +- .../gonum.org/v1/gonum/lapack/gonum/dorgqr.go | 8 +- .../gonum.org/v1/gonum/lapack/gonum/dorgtr.go | 2 +- .../gonum.org/v1/gonum/lapack/gonum/dorm2r.go | 6 +- .../gonum.org/v1/gonum/lapack/gonum/dptcon.go | 99 +++ .../gonum.org/v1/gonum/lapack/gonum/dptsv.go | 49 ++ .../gonum.org/v1/gonum/lapack/gonum/dpttrf.go | 80 ++ .../gonum.org/v1/gonum/lapack/gonum/dpttrs.go | 51 ++ .../gonum.org/v1/gonum/lapack/gonum/dptts2.go | 39 + .../gonum.org/v1/gonum/lapack/gonum/errors.go | 1 + .../gonum.org/v1/gonum/lapack/gonum/ilaenv.go | 7 + .../gonum.org/v1/gonum/lapack/gonum/lapack.go | 14 - vendor/gonum.org/v1/gonum/lapack/lapack.go | 14 +- .../v1/gonum/lapack/lapack64/lapack64.go | 146 +++- vendor/gonum.org/v1/gonum/mat/cholesky.go | 309 +++++++- vendor/gonum.org/v1/gonum/mat/dense.go | 55 ++ vendor/gonum.org/v1/gonum/mat/eigen.go | 125 +++- vendor/gonum.org/v1/gonum/mat/lq.go | 80 +- vendor/gonum.org/v1/gonum/mat/lu.go | 242 +++--- vendor/gonum.org/v1/gonum/mat/matrix.go | 25 +- vendor/gonum.org/v1/gonum/mat/qr.go | 108 ++- vendor/gonum.org/v1/gonum/mat/solve.go | 62 +- vendor/gonum.org/v1/gonum/mat/triangular.go | 51 ++ vendor/gonum.org/v1/gonum/mat/triband.go | 18 +- vendor/gonum.org/v1/gonum/mat/tridiag.go | 18 +- vendor/gonum.org/v1/gonum/mat/vector.go | 16 + vendor/modules.txt | 6 +- 66 files changed, 2195 insertions(+), 1049 deletions(-) rename vendor/github.com/google/go-cmp/cmp/{export_unsafe.go => export.go} (94%) delete mode 100644 vendor/github.com/google/go-cmp/cmp/export_panic.go rename vendor/github.com/google/go-cmp/cmp/internal/value/{pointer_unsafe.go => pointer.go} (95%) delete mode 100644 vendor/github.com/google/go-cmp/cmp/internal/value/pointer_purego.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dgghrd.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dlanhs.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dptcon.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dptsv.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dpttrf.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dpttrs.go create mode 100644 vendor/gonum.org/v1/gonum/lapack/gonum/dptts2.go diff --git a/go.mod b/go.mod index 4eacd5c9..434c3f64 100644 --- a/go.mod +++ b/go.mod @@ -13,7 +13,7 @@ require ( github.com/facebookgo/grace v0.0.0-20180706040059-75cf19382434 github.com/facebookgo/pidfile v0.0.0-20150612191647-f242e2999868 github.com/go-graphite/protocol v1.0.1-0.20220718132526-4b842ba389ee - github.com/google/go-cmp v0.5.9 + github.com/google/go-cmp v0.6.0 github.com/gorilla/handlers v1.5.1 github.com/gorilla/mux v1.8.0 github.com/lomik/og-rek v0.0.0-20170411191824-628eefeb8d80 @@ -25,7 +25,7 @@ require ( github.com/wangjohn/quickselect v0.0.0-20161129230411-ed8402a42d5f go.opentelemetry.io/contrib/instrumentation/gorilla/mux v0.7.0 go.uber.org/zap v1.24.0 - gonum.org/v1/gonum v0.13.0 + gonum.org/v1/gonum v0.15.1 google.golang.org/grpc v1.56.3 gopkg.in/yaml.v2 v2.4.0 ) diff --git a/go.sum b/go.sum index 8dd8a4d5..182b4adb 100644 --- a/go.sum +++ b/go.sum @@ -93,8 +93,8 @@ github.com/google/go-cmp v0.4.0/go.mod h1:v8dTdLbMG2kIc/vJvl+f65V22dbkXbowE6jgT/ github.com/google/go-cmp v0.5.0/go.mod h1:v8dTdLbMG2kIc/vJvl+f65V22dbkXbowE6jgT/gNBxE= github.com/google/go-cmp v0.5.5/go.mod h1:v8dTdLbMG2kIc/vJvl+f65V22dbkXbowE6jgT/gNBxE= github.com/google/go-cmp v0.5.6/go.mod h1:v8dTdLbMG2kIc/vJvl+f65V22dbkXbowE6jgT/gNBxE= -github.com/google/go-cmp v0.5.9 h1:O2Tfq5qg4qc4AmwVlvv0oLiVAGB7enBSJ2x2DqQFi38= -github.com/google/go-cmp v0.5.9/go.mod h1:17dUlkBOakJ0+DkrSSNjCkIjxS6bF9zb3elmeNGIjoY= +github.com/google/go-cmp v0.6.0 h1:ofyhxvXcZhMsU5ulbFiLKl/XBFqE1GSq7atu8tAmTRI= +github.com/google/go-cmp v0.6.0/go.mod h1:17dUlkBOakJ0+DkrSSNjCkIjxS6bF9zb3elmeNGIjoY= github.com/google/gofuzz v1.0.0/go.mod h1:dBl0BpW6vV/+mYPU4Po3pmUjxk6FQPldtuIdl/M65Eg= github.com/google/uuid v1.1.2/go.mod h1:TIyPZe4MgqvfeYDBFedMoGGpEw/LqOeaOT+nhxU+yHo= github.com/gorilla/handlers v1.5.1 h1:9lRY6j8DEeeBT10CvO9hGW0gmky0BprnvDI5vfhUHH4= @@ -162,7 +162,7 @@ go.uber.org/zap v1.24.0/go.mod h1:2kMP+WWQ8aoFoedH3T2sq6iJ2yDWpHbP0f6MQbS9Gkg= golang.org/x/crypto v0.0.0-20190308221718-c2843e01d9a2/go.mod h1:djNgcEr1/C05ACkg1iLfiJU5Ep61QUkGW8qpdssI0+w= golang.org/x/crypto v0.0.0-20200622213623-75b288015ac9/go.mod h1:LzIPMQfyMNhhGPhUkYOs5KpL4U8rLKemX1yGLhDgUto= golang.org/x/exp v0.0.0-20190121172915-509febef88a4/go.mod h1:CJ0aWSM057203Lf6IL+f9T1iT9GByDxfZKAQTCR3kQA= -golang.org/x/exp v0.0.0-20230321023759-10a507213a29 h1:ooxPy7fPvB4kwsA2h+iBNHkAbp/4JxTSwCmvdjEYmug= +golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa h1:FRnLl4eNAQl8hwxVVC17teOw8kdjVDVAiFMtgUdTSRQ= golang.org/x/lint v0.0.0-20181026193005-c67002cb31c3/go.mod h1:UVdnD1Gm6xHRNCYTkRU2/jEulfH38KcIWyp/GAMgvoE= golang.org/x/lint v0.0.0-20190227174305-5b3e6a55c961/go.mod h1:wehouNa3lNwaWXcvxsM5YxQ5yQlVC4a0KAMCusXpPoU= golang.org/x/lint v0.0.0-20190313153728-d0100b6bd8b3/go.mod h1:6SW0HCj/g11FgYtHlgUYUwCkIfeOF89ocIRzGO/8vkc= @@ -203,8 +203,8 @@ golang.org/x/tools v0.0.0-20190311212946-11955173bddd/go.mod h1:LCzVGOaR6xXOjkQ3 golang.org/x/tools v0.0.0-20190524140312-2c0ae7006135/go.mod h1:RgjU9mgBXZiqYHBnxXauZ1Gv1EHHAz9KjViQ78xBX0Q= golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= golang.org/x/xerrors v0.0.0-20200804184101-5ec99f83aff1/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= -gonum.org/v1/gonum v0.13.0 h1:a0T3bh+7fhRyqeNbiC3qVHYmkiQgit3wnNan/2c0HMM= -gonum.org/v1/gonum v0.13.0/go.mod h1:/WPYRckkfWrhWefxyYTfrTtQR0KH4iyHNuzxqXAKyAU= +gonum.org/v1/gonum v0.15.1 h1:FNy7N6OUZVUaWG9pTiD+jlhdQ3lMP+/LcTpJ6+a8sQ0= +gonum.org/v1/gonum v0.15.1/go.mod h1:eZTZuRFrzu5pcyjN5wJhcIhnUdNijYxX1T2IcrOGY0o= google.golang.org/appengine v1.1.0/go.mod h1:EbEs0AVv82hx2wNQdGPgUI5lhzA/G0D9YwlJXL52JkM= google.golang.org/appengine v1.4.0/go.mod h1:xpcJRLb0r/rnEns0DIKYYv+WjYCduHsrkT7/EB5XEv4= google.golang.org/genproto v0.0.0-20180817151627-c66870c02cf8/go.mod h1:JiN7NxoALGmiZfu7CAH4rXhgtRTLTxftemlI0sWmxmc= diff --git a/vendor/github.com/google/go-cmp/cmp/compare.go b/vendor/github.com/google/go-cmp/cmp/compare.go index 087320da..0f5b8a48 100644 --- a/vendor/github.com/google/go-cmp/cmp/compare.go +++ b/vendor/github.com/google/go-cmp/cmp/compare.go @@ -5,7 +5,7 @@ // Package cmp determines equality of values. // // This package is intended to be a more powerful and safer alternative to -// reflect.DeepEqual for comparing whether two values are semantically equal. +// [reflect.DeepEqual] for comparing whether two values are semantically equal. // It is intended to only be used in tests, as performance is not a goal and // it may panic if it cannot compare the values. Its propensity towards // panicking means that its unsuitable for production environments where a @@ -18,16 +18,17 @@ // For example, an equality function may report floats as equal so long as // they are within some tolerance of each other. // -// - Types with an Equal method may use that method to determine equality. -// This allows package authors to determine the equality operation -// for the types that they define. +// - Types with an Equal method (e.g., [time.Time.Equal]) may use that method +// to determine equality. This allows package authors to determine +// the equality operation for the types that they define. // // - If no custom equality functions are used and no Equal method is defined, // equality is determined by recursively comparing the primitive kinds on -// both values, much like reflect.DeepEqual. Unlike reflect.DeepEqual, +// both values, much like [reflect.DeepEqual]. Unlike [reflect.DeepEqual], // unexported fields are not compared by default; they result in panics -// unless suppressed by using an Ignore option (see cmpopts.IgnoreUnexported) -// or explicitly compared using the Exporter option. +// unless suppressed by using an [Ignore] option +// (see [github.com/google/go-cmp/cmp/cmpopts.IgnoreUnexported]) +// or explicitly compared using the [Exporter] option. package cmp import ( @@ -45,14 +46,14 @@ import ( // Equal reports whether x and y are equal by recursively applying the // following rules in the given order to x and y and all of their sub-values: // -// - Let S be the set of all Ignore, Transformer, and Comparer options that +// - Let S be the set of all [Ignore], [Transformer], and [Comparer] options that // remain after applying all path filters, value filters, and type filters. -// If at least one Ignore exists in S, then the comparison is ignored. -// If the number of Transformer and Comparer options in S is non-zero, +// If at least one [Ignore] exists in S, then the comparison is ignored. +// If the number of [Transformer] and [Comparer] options in S is non-zero, // then Equal panics because it is ambiguous which option to use. -// If S contains a single Transformer, then use that to transform +// If S contains a single [Transformer], then use that to transform // the current values and recursively call Equal on the output values. -// If S contains a single Comparer, then use that to compare the current values. +// If S contains a single [Comparer], then use that to compare the current values. // Otherwise, evaluation proceeds to the next rule. // // - If the values have an Equal method of the form "(T) Equal(T) bool" or @@ -66,21 +67,22 @@ import ( // Functions are only equal if they are both nil, otherwise they are unequal. // // Structs are equal if recursively calling Equal on all fields report equal. -// If a struct contains unexported fields, Equal panics unless an Ignore option -// (e.g., cmpopts.IgnoreUnexported) ignores that field or the Exporter option -// explicitly permits comparing the unexported field. +// If a struct contains unexported fields, Equal panics unless an [Ignore] option +// (e.g., [github.com/google/go-cmp/cmp/cmpopts.IgnoreUnexported]) ignores that field +// or the [Exporter] option explicitly permits comparing the unexported field. // // Slices are equal if they are both nil or both non-nil, where recursively // calling Equal on all non-ignored slice or array elements report equal. // Empty non-nil slices and nil slices are not equal; to equate empty slices, -// consider using cmpopts.EquateEmpty. +// consider using [github.com/google/go-cmp/cmp/cmpopts.EquateEmpty]. // // Maps are equal if they are both nil or both non-nil, where recursively // calling Equal on all non-ignored map entries report equal. // Map keys are equal according to the == operator. -// To use custom comparisons for map keys, consider using cmpopts.SortMaps. +// To use custom comparisons for map keys, consider using +// [github.com/google/go-cmp/cmp/cmpopts.SortMaps]. // Empty non-nil maps and nil maps are not equal; to equate empty maps, -// consider using cmpopts.EquateEmpty. +// consider using [github.com/google/go-cmp/cmp/cmpopts.EquateEmpty]. // // Pointers and interfaces are equal if they are both nil or both non-nil, // where they have the same underlying concrete type and recursively diff --git a/vendor/github.com/google/go-cmp/cmp/export_unsafe.go b/vendor/github.com/google/go-cmp/cmp/export.go similarity index 94% rename from vendor/github.com/google/go-cmp/cmp/export_unsafe.go rename to vendor/github.com/google/go-cmp/cmp/export.go index e2c0f74e..29f82fe6 100644 --- a/vendor/github.com/google/go-cmp/cmp/export_unsafe.go +++ b/vendor/github.com/google/go-cmp/cmp/export.go @@ -2,9 +2,6 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -//go:build !purego -// +build !purego - package cmp import ( @@ -12,8 +9,6 @@ import ( "unsafe" ) -const supportExporters = true - // retrieveUnexportedField uses unsafe to forcibly retrieve any field from // a struct such that the value has read-write permissions. // diff --git a/vendor/github.com/google/go-cmp/cmp/export_panic.go b/vendor/github.com/google/go-cmp/cmp/export_panic.go deleted file mode 100644 index ae851fe5..00000000 --- a/vendor/github.com/google/go-cmp/cmp/export_panic.go +++ /dev/null @@ -1,16 +0,0 @@ -// Copyright 2017, The Go Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -//go:build purego -// +build purego - -package cmp - -import "reflect" - -const supportExporters = false - -func retrieveUnexportedField(reflect.Value, reflect.StructField, bool) reflect.Value { - panic("no support for forcibly accessing unexported fields") -} diff --git a/vendor/github.com/google/go-cmp/cmp/internal/value/pointer_unsafe.go b/vendor/github.com/google/go-cmp/cmp/internal/value/pointer.go similarity index 95% rename from vendor/github.com/google/go-cmp/cmp/internal/value/pointer_unsafe.go rename to vendor/github.com/google/go-cmp/cmp/internal/value/pointer.go index 16e6860a..e5dfff69 100644 --- a/vendor/github.com/google/go-cmp/cmp/internal/value/pointer_unsafe.go +++ b/vendor/github.com/google/go-cmp/cmp/internal/value/pointer.go @@ -2,9 +2,6 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -//go:build !purego -// +build !purego - package value import ( diff --git a/vendor/github.com/google/go-cmp/cmp/internal/value/pointer_purego.go b/vendor/github.com/google/go-cmp/cmp/internal/value/pointer_purego.go deleted file mode 100644 index 1a71bfcb..00000000 --- a/vendor/github.com/google/go-cmp/cmp/internal/value/pointer_purego.go +++ /dev/null @@ -1,34 +0,0 @@ -// Copyright 2018, The Go Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -//go:build purego -// +build purego - -package value - -import "reflect" - -// Pointer is an opaque typed pointer and is guaranteed to be comparable. -type Pointer struct { - p uintptr - t reflect.Type -} - -// PointerOf returns a Pointer from v, which must be a -// reflect.Ptr, reflect.Slice, or reflect.Map. -func PointerOf(v reflect.Value) Pointer { - // NOTE: Storing a pointer as an uintptr is technically incorrect as it - // assumes that the GC implementation does not use a moving collector. - return Pointer{v.Pointer(), v.Type()} -} - -// IsNil reports whether the pointer is nil. -func (p Pointer) IsNil() bool { - return p.p == 0 -} - -// Uintptr returns the pointer as a uintptr. -func (p Pointer) Uintptr() uintptr { - return p.p -} diff --git a/vendor/github.com/google/go-cmp/cmp/options.go b/vendor/github.com/google/go-cmp/cmp/options.go index 1f9ca9c4..754496f3 100644 --- a/vendor/github.com/google/go-cmp/cmp/options.go +++ b/vendor/github.com/google/go-cmp/cmp/options.go @@ -13,15 +13,15 @@ import ( "github.com/google/go-cmp/cmp/internal/function" ) -// Option configures for specific behavior of Equal and Diff. In particular, -// the fundamental Option functions (Ignore, Transformer, and Comparer), +// Option configures for specific behavior of [Equal] and [Diff]. In particular, +// the fundamental Option functions ([Ignore], [Transformer], and [Comparer]), // configure how equality is determined. // -// The fundamental options may be composed with filters (FilterPath and -// FilterValues) to control the scope over which they are applied. +// The fundamental options may be composed with filters ([FilterPath] and +// [FilterValues]) to control the scope over which they are applied. // -// The cmp/cmpopts package provides helper functions for creating options that -// may be used with Equal and Diff. +// The [github.com/google/go-cmp/cmp/cmpopts] package provides helper functions +// for creating options that may be used with [Equal] and [Diff]. type Option interface { // filter applies all filters and returns the option that remains. // Each option may only read s.curPath and call s.callTTBFunc. @@ -56,9 +56,9 @@ type core struct{} func (core) isCore() {} -// Options is a list of Option values that also satisfies the Option interface. +// Options is a list of [Option] values that also satisfies the [Option] interface. // Helper comparison packages may return an Options value when packing multiple -// Option values into a single Option. When this package processes an Options, +// [Option] values into a single [Option]. When this package processes an Options, // it will be implicitly expanded into a flat list. // // Applying a filter on an Options is equivalent to applying that same filter @@ -105,16 +105,16 @@ func (opts Options) String() string { return fmt.Sprintf("Options{%s}", strings.Join(ss, ", ")) } -// FilterPath returns a new Option where opt is only evaluated if filter f -// returns true for the current Path in the value tree. +// FilterPath returns a new [Option] where opt is only evaluated if filter f +// returns true for the current [Path] in the value tree. // // This filter is called even if a slice element or map entry is missing and // provides an opportunity to ignore such cases. The filter function must be // symmetric such that the filter result is identical regardless of whether the // missing value is from x or y. // -// The option passed in may be an Ignore, Transformer, Comparer, Options, or -// a previously filtered Option. +// The option passed in may be an [Ignore], [Transformer], [Comparer], [Options], or +// a previously filtered [Option]. func FilterPath(f func(Path) bool, opt Option) Option { if f == nil { panic("invalid path filter function") @@ -142,7 +142,7 @@ func (f pathFilter) String() string { return fmt.Sprintf("FilterPath(%s, %v)", function.NameOf(reflect.ValueOf(f.fnc)), f.opt) } -// FilterValues returns a new Option where opt is only evaluated if filter f, +// FilterValues returns a new [Option] where opt is only evaluated if filter f, // which is a function of the form "func(T, T) bool", returns true for the // current pair of values being compared. If either value is invalid or // the type of the values is not assignable to T, then this filter implicitly @@ -154,8 +154,8 @@ func (f pathFilter) String() string { // If T is an interface, it is possible that f is called with two values with // different concrete types that both implement T. // -// The option passed in may be an Ignore, Transformer, Comparer, Options, or -// a previously filtered Option. +// The option passed in may be an [Ignore], [Transformer], [Comparer], [Options], or +// a previously filtered [Option]. func FilterValues(f interface{}, opt Option) Option { v := reflect.ValueOf(f) if !function.IsType(v.Type(), function.ValueFilter) || v.IsNil() { @@ -192,9 +192,9 @@ func (f valuesFilter) String() string { return fmt.Sprintf("FilterValues(%s, %v)", function.NameOf(f.fnc), f.opt) } -// Ignore is an Option that causes all comparisons to be ignored. -// This value is intended to be combined with FilterPath or FilterValues. -// It is an error to pass an unfiltered Ignore option to Equal. +// Ignore is an [Option] that causes all comparisons to be ignored. +// This value is intended to be combined with [FilterPath] or [FilterValues]. +// It is an error to pass an unfiltered Ignore option to [Equal]. func Ignore() Option { return ignore{} } type ignore struct{ core } @@ -234,6 +234,8 @@ func (validator) apply(s *state, vx, vy reflect.Value) { name = fmt.Sprintf("%q.%v", t.PkgPath(), t.Name()) // e.g., "path/to/package".MyType if _, ok := reflect.New(t).Interface().(error); ok { help = "consider using cmpopts.EquateErrors to compare error values" + } else if t.Comparable() { + help = "consider using cmpopts.EquateComparable to compare comparable Go types" } } else { // Unnamed type with unexported fields. Derive PkgPath from field. @@ -254,7 +256,7 @@ const identRx = `[_\p{L}][_\p{L}\p{N}]*` var identsRx = regexp.MustCompile(`^` + identRx + `(\.` + identRx + `)*$`) -// Transformer returns an Option that applies a transformation function that +// Transformer returns an [Option] that applies a transformation function that // converts values of a certain type into that of another. // // The transformer f must be a function "func(T) R" that converts values of @@ -265,13 +267,14 @@ var identsRx = regexp.MustCompile(`^` + identRx + `(\.` + identRx + `)*$`) // same transform to the output of itself (e.g., in the case where the // input and output types are the same), an implicit filter is added such that // a transformer is applicable only if that exact transformer is not already -// in the tail of the Path since the last non-Transform step. +// in the tail of the [Path] since the last non-[Transform] step. // For situations where the implicit filter is still insufficient, -// consider using cmpopts.AcyclicTransformer, which adds a filter -// to prevent the transformer from being recursively applied upon itself. +// consider using [github.com/google/go-cmp/cmp/cmpopts.AcyclicTransformer], +// which adds a filter to prevent the transformer from +// being recursively applied upon itself. // -// The name is a user provided label that is used as the Transform.Name in the -// transformation PathStep (and eventually shown in the Diff output). +// The name is a user provided label that is used as the [Transform.Name] in the +// transformation [PathStep] (and eventually shown in the [Diff] output). // The name must be a valid identifier or qualified identifier in Go syntax. // If empty, an arbitrary name is used. func Transformer(name string, f interface{}) Option { @@ -329,7 +332,7 @@ func (tr transformer) String() string { return fmt.Sprintf("Transformer(%s, %s)", tr.name, function.NameOf(tr.fnc)) } -// Comparer returns an Option that determines whether two values are equal +// Comparer returns an [Option] that determines whether two values are equal // to each other. // // The comparer f must be a function "func(T, T) bool" and is implicitly @@ -377,35 +380,32 @@ func (cm comparer) String() string { return fmt.Sprintf("Comparer(%s)", function.NameOf(cm.fnc)) } -// Exporter returns an Option that specifies whether Equal is allowed to +// Exporter returns an [Option] that specifies whether [Equal] is allowed to // introspect into the unexported fields of certain struct types. // // Users of this option must understand that comparing on unexported fields // from external packages is not safe since changes in the internal -// implementation of some external package may cause the result of Equal +// implementation of some external package may cause the result of [Equal] // to unexpectedly change. However, it may be valid to use this option on types // defined in an internal package where the semantic meaning of an unexported // field is in the control of the user. // -// In many cases, a custom Comparer should be used instead that defines +// In many cases, a custom [Comparer] should be used instead that defines // equality as a function of the public API of a type rather than the underlying // unexported implementation. // -// For example, the reflect.Type documentation defines equality to be determined +// For example, the [reflect.Type] documentation defines equality to be determined // by the == operator on the interface (essentially performing a shallow pointer -// comparison) and most attempts to compare *regexp.Regexp types are interested +// comparison) and most attempts to compare *[regexp.Regexp] types are interested // in only checking that the regular expression strings are equal. -// Both of these are accomplished using Comparers: +// Both of these are accomplished using [Comparer] options: // // Comparer(func(x, y reflect.Type) bool { return x == y }) // Comparer(func(x, y *regexp.Regexp) bool { return x.String() == y.String() }) // -// In other cases, the cmpopts.IgnoreUnexported option can be used to ignore -// all unexported fields on specified struct types. +// In other cases, the [github.com/google/go-cmp/cmp/cmpopts.IgnoreUnexported] +// option can be used to ignore all unexported fields on specified struct types. func Exporter(f func(reflect.Type) bool) Option { - if !supportExporters { - panic("Exporter is not supported on purego builds") - } return exporter(f) } @@ -415,10 +415,10 @@ func (exporter) filter(_ *state, _ reflect.Type, _, _ reflect.Value) applicableO panic("not implemented") } -// AllowUnexported returns an Options that allows Equal to forcibly introspect +// AllowUnexported returns an [Option] that allows [Equal] to forcibly introspect // unexported fields of the specified struct types. // -// See Exporter for the proper use of this option. +// See [Exporter] for the proper use of this option. func AllowUnexported(types ...interface{}) Option { m := make(map[reflect.Type]bool) for _, typ := range types { @@ -432,7 +432,7 @@ func AllowUnexported(types ...interface{}) Option { } // Result represents the comparison result for a single node and -// is provided by cmp when calling Report (see Reporter). +// is provided by cmp when calling Report (see [Reporter]). type Result struct { _ [0]func() // Make Result incomparable flags resultFlags @@ -445,7 +445,7 @@ func (r Result) Equal() bool { } // ByIgnore reports whether the node is equal because it was ignored. -// This never reports true if Equal reports false. +// This never reports true if [Result.Equal] reports false. func (r Result) ByIgnore() bool { return r.flags&reportByIgnore != 0 } @@ -455,7 +455,7 @@ func (r Result) ByMethod() bool { return r.flags&reportByMethod != 0 } -// ByFunc reports whether a Comparer function determined equality. +// ByFunc reports whether a [Comparer] function determined equality. func (r Result) ByFunc() bool { return r.flags&reportByFunc != 0 } @@ -478,7 +478,7 @@ const ( reportByCycle ) -// Reporter is an Option that can be passed to Equal. When Equal traverses +// Reporter is an [Option] that can be passed to [Equal]. When [Equal] traverses // the value trees, it calls PushStep as it descends into each node in the // tree and PopStep as it ascend out of the node. The leaves of the tree are // either compared (determined to be equal or not equal) or ignored and reported diff --git a/vendor/github.com/google/go-cmp/cmp/path.go b/vendor/github.com/google/go-cmp/cmp/path.go index a0a58850..c3c14564 100644 --- a/vendor/github.com/google/go-cmp/cmp/path.go +++ b/vendor/github.com/google/go-cmp/cmp/path.go @@ -14,9 +14,9 @@ import ( "github.com/google/go-cmp/cmp/internal/value" ) -// Path is a list of PathSteps describing the sequence of operations to get +// Path is a list of [PathStep] describing the sequence of operations to get // from some root type to the current position in the value tree. -// The first Path element is always an operation-less PathStep that exists +// The first Path element is always an operation-less [PathStep] that exists // simply to identify the initial type. // // When traversing structs with embedded structs, the embedded struct will @@ -29,8 +29,13 @@ type Path []PathStep // a value's tree structure. Users of this package never need to implement // these types as values of this type will be returned by this package. // -// Implementations of this interface are -// StructField, SliceIndex, MapIndex, Indirect, TypeAssertion, and Transform. +// Implementations of this interface: +// - [StructField] +// - [SliceIndex] +// - [MapIndex] +// - [Indirect] +// - [TypeAssertion] +// - [Transform] type PathStep interface { String() string @@ -70,8 +75,9 @@ func (pa *Path) pop() { *pa = (*pa)[:len(*pa)-1] } -// Last returns the last PathStep in the Path. -// If the path is empty, this returns a non-nil PathStep that reports a nil Type. +// Last returns the last [PathStep] in the Path. +// If the path is empty, this returns a non-nil [PathStep] +// that reports a nil [PathStep.Type]. func (pa Path) Last() PathStep { return pa.Index(-1) } @@ -79,7 +85,8 @@ func (pa Path) Last() PathStep { // Index returns the ith step in the Path and supports negative indexing. // A negative index starts counting from the tail of the Path such that -1 // refers to the last step, -2 refers to the second-to-last step, and so on. -// If index is invalid, this returns a non-nil PathStep that reports a nil Type. +// If index is invalid, this returns a non-nil [PathStep] +// that reports a nil [PathStep.Type]. func (pa Path) Index(i int) PathStep { if i < 0 { i = len(pa) + i @@ -168,7 +175,8 @@ func (ps pathStep) String() string { return fmt.Sprintf("{%s}", s) } -// StructField represents a struct field access on a field called Name. +// StructField is a [PathStep] that represents a struct field access +// on a field called [StructField.Name]. type StructField struct{ *structField } type structField struct { pathStep @@ -204,10 +212,11 @@ func (sf StructField) String() string { return fmt.Sprintf(".%s", sf.name) } func (sf StructField) Name() string { return sf.name } // Index is the index of the field in the parent struct type. -// See reflect.Type.Field. +// See [reflect.Type.Field]. func (sf StructField) Index() int { return sf.idx } -// SliceIndex is an index operation on a slice or array at some index Key. +// SliceIndex is a [PathStep] that represents an index operation on +// a slice or array at some index [SliceIndex.Key]. type SliceIndex struct{ *sliceIndex } type sliceIndex struct { pathStep @@ -247,12 +256,12 @@ func (si SliceIndex) Key() int { // all of the indexes to be shifted. If an index is -1, then that // indicates that the element does not exist in the associated slice. // -// Key is guaranteed to return -1 if and only if the indexes returned -// by SplitKeys are not the same. SplitKeys will never return -1 for +// [SliceIndex.Key] is guaranteed to return -1 if and only if the indexes +// returned by SplitKeys are not the same. SplitKeys will never return -1 for // both indexes. func (si SliceIndex) SplitKeys() (ix, iy int) { return si.xkey, si.ykey } -// MapIndex is an index operation on a map at some index Key. +// MapIndex is a [PathStep] that represents an index operation on a map at some index Key. type MapIndex struct{ *mapIndex } type mapIndex struct { pathStep @@ -266,7 +275,7 @@ func (mi MapIndex) String() string { return fmt.Sprintf("[%#v]", // Key is the value of the map key. func (mi MapIndex) Key() reflect.Value { return mi.key } -// Indirect represents pointer indirection on the parent type. +// Indirect is a [PathStep] that represents pointer indirection on the parent type. type Indirect struct{ *indirect } type indirect struct { pathStep @@ -276,7 +285,7 @@ func (in Indirect) Type() reflect.Type { return in.typ } func (in Indirect) Values() (vx, vy reflect.Value) { return in.vx, in.vy } func (in Indirect) String() string { return "*" } -// TypeAssertion represents a type assertion on an interface. +// TypeAssertion is a [PathStep] that represents a type assertion on an interface. type TypeAssertion struct{ *typeAssertion } type typeAssertion struct { pathStep @@ -286,7 +295,8 @@ func (ta TypeAssertion) Type() reflect.Type { return ta.typ } func (ta TypeAssertion) Values() (vx, vy reflect.Value) { return ta.vx, ta.vy } func (ta TypeAssertion) String() string { return fmt.Sprintf(".(%v)", value.TypeString(ta.typ, false)) } -// Transform is a transformation from the parent type to the current type. +// Transform is a [PathStep] that represents a transformation +// from the parent type to the current type. type Transform struct{ *transform } type transform struct { pathStep @@ -297,13 +307,13 @@ func (tf Transform) Type() reflect.Type { return tf.typ } func (tf Transform) Values() (vx, vy reflect.Value) { return tf.vx, tf.vy } func (tf Transform) String() string { return fmt.Sprintf("%s()", tf.trans.name) } -// Name is the name of the Transformer. +// Name is the name of the [Transformer]. func (tf Transform) Name() string { return tf.trans.name } // Func is the function pointer to the transformer function. func (tf Transform) Func() reflect.Value { return tf.trans.fnc } -// Option returns the originally constructed Transformer option. +// Option returns the originally constructed [Transformer] option. // The == operator can be used to detect the exact option used. func (tf Transform) Option() Option { return tf.trans } diff --git a/vendor/github.com/google/go-cmp/cmp/report_reflect.go b/vendor/github.com/google/go-cmp/cmp/report_reflect.go index 2ab41fad..e39f4228 100644 --- a/vendor/github.com/google/go-cmp/cmp/report_reflect.go +++ b/vendor/github.com/google/go-cmp/cmp/report_reflect.go @@ -199,7 +199,7 @@ func (opts formatOptions) FormatValue(v reflect.Value, parentKind reflect.Kind, break } sf := t.Field(i) - if supportExporters && !isExported(sf.Name) { + if !isExported(sf.Name) { vv = retrieveUnexportedField(v, sf, true) } s := opts.WithTypeMode(autoType).FormatValue(vv, t.Kind(), ptrs) diff --git a/vendor/gonum.org/v1/gonum/AUTHORS b/vendor/gonum.org/v1/gonum/AUTHORS index 7d49714a..8e5896db 100644 --- a/vendor/gonum.org/v1/gonum/AUTHORS +++ b/vendor/gonum.org/v1/gonum/AUTHORS @@ -11,10 +11,12 @@ Alexander Egurnov Andrei Blinnikov antichris +Bailey Lissington Bill Gray Bill Noon Brendan Tracey Brent Pedersen +Bulat Khasanov Chad Kunde Chan Kwan Yin Chih-Wei Chang @@ -32,11 +34,14 @@ Davor Kapsa DeepMind Technologies Delaney Gillilan Dezmond Goff +Dirk Müller Dong-hee Na Dustin Spicuzza Egon Elbre Ekaterina Efimova +Eng Zer Jun Ethan Burns +Ethan Reesor Evert Lammerts Evgeny Savinov Fabian Wickborn @@ -46,6 +51,7 @@ Francesc Campoy Google Inc Gustaf Johansson Hossein Zolfi +Huang Peng Fei Iakov Davydov Igor Mikushkin Iskander Sharipov @@ -55,17 +61,21 @@ James Bowman James Holmes <32bitkid@gmail.com> Janne Snabb Jeremy Atkinson +Jes Cok Jinesi Yelizati Jonas Kahler Jonas Schulze +Jonathan Bluett-Duncan Jonathan J Lawlor Jonathan Reiter Jonathan Schroeder Joost van Amersfoort +Jordan Stoker Joseph Watson Josh Wilson Julien Roland Kai Trukenmüller +Kendall Marcus Kent English Kevin C. Zimmerman Kirill Motkov @@ -114,6 +124,8 @@ The University of Minnesota The University of Washington Thomas Berg Tobin Harding +Tom Payne +Tristan Nicholls Valentin Deleplace Vincent Thiery Vladimír Chalupecký diff --git a/vendor/gonum.org/v1/gonum/CONTRIBUTORS b/vendor/gonum.org/v1/gonum/CONTRIBUTORS index b8bef3e3..e367595b 100644 --- a/vendor/gonum.org/v1/gonum/CONTRIBUTORS +++ b/vendor/gonum.org/v1/gonum/CONTRIBUTORS @@ -19,10 +19,12 @@ Alexander Egurnov Andrei Blinnikov Andrew Brampton antichris +Bailey Lissington Bill Gray Bill Noon Brendan Tracey Brent Pedersen +Bulat Khasanov Chad Kunde Chan Kwan Yin Chih-Wei Chang @@ -40,11 +42,14 @@ David Samborski Davor Kapsa Delaney Gillilan Dezmond Goff +Dirk Müller Dong-hee Na Dustin Spicuzza Egon Elbre Ekaterina Efimova +Eng Zer Jun Ethan Burns +Ethan Reesor Evert Lammerts Evgeny Savinov Fabian Wickborn @@ -53,6 +58,7 @@ Fazlul Shahriar Francesc Campoy Gustaf Johansson Hossein Zolfi +Huang Peng Fei Iakov Davydov Igor Mikushkin Iskander Sharipov @@ -62,18 +68,22 @@ James Bowman James Holmes <32bitkid@gmail.com> Janne Snabb Jeremy Atkinson +Jes Cok Jinesi Yelizati Jon Richards Jonas Kahler Jonas Schulze +Jonathan Bluett-Duncan Jonathan J Lawlor Jonathan Reiter Jonathan Schroeder Joost van Amersfoort +Jordan Stoker Joseph Watson Josh Wilson Julien Roland Kai Trukenmüller +Kendall Marcus Kent English Kevin C. Zimmerman Kirill Motkov @@ -117,6 +127,8 @@ Takeshi Yoneda Tamir Hyman Thomas Berg Tobin Harding +Tom Payne +Tristan Nicholls Valentin Deleplace Vincent Thiery Vladimír Chalupecký diff --git a/vendor/gonum.org/v1/gonum/blas/blas64/blas64.go b/vendor/gonum.org/v1/gonum/blas/blas64/blas64.go index c336dc89..64ac985c 100644 --- a/vendor/gonum.org/v1/gonum/blas/blas64/blas64.go +++ b/vendor/gonum.org/v1/gonum/blas/blas64/blas64.go @@ -20,7 +20,7 @@ func Use(b blas.Float64) { // Implementation returns the current BLAS float64 implementation. // -// Implementation allows direct calls to the current the BLAS float64 implementation +// Implementation allows direct calls to the current BLAS float64 implementation // giving finer control of parameters. func Implementation() blas.Float64 { return blas64 diff --git a/vendor/gonum.org/v1/gonum/blas/blas64/conv.go b/vendor/gonum.org/v1/gonum/blas/blas64/conv.go index 6cc6517f..695557d1 100644 --- a/vendor/gonum.org/v1/gonum/blas/blas64/conv.go +++ b/vendor/gonum.org/v1/gonum/blas/blas64/conv.go @@ -261,17 +261,3 @@ func (t TriangularBand) From(a TriangularBandCols) { } dst.From(src) } - -func min(a, b int) int { - if a < b { - return a - } - return b -} - -func max(a, b int) int { - if a > b { - return a - } - return b -} diff --git a/vendor/gonum.org/v1/gonum/blas/cblas128/conv.go b/vendor/gonum.org/v1/gonum/blas/cblas128/conv.go index c459e1d8..bfafb96e 100644 --- a/vendor/gonum.org/v1/gonum/blas/cblas128/conv.go +++ b/vendor/gonum.org/v1/gonum/blas/cblas128/conv.go @@ -263,17 +263,3 @@ func (t TriangularBand) From(a TriangularBandCols) { } dst.From(src) } - -func min(a, b int) int { - if a < b { - return a - } - return b -} - -func max(a, b int) int { - if a > b { - return a - } - return b -} diff --git a/vendor/gonum.org/v1/gonum/blas/gonum/gonum.go b/vendor/gonum.org/v1/gonum/blas/gonum/gonum.go index 61a8b8b5..5a5c1110 100644 --- a/vendor/gonum.org/v1/gonum/blas/gonum/gonum.go +++ b/vendor/gonum.org/v1/gonum/blas/gonum/gonum.go @@ -21,20 +21,6 @@ const ( minParBlock = 4 // minimum number of blocks needed to go parallel ) -func max(a, b int) int { - if a > b { - return a - } - return b -} - -func min(a, b int) int { - if a > b { - return b - } - return a -} - // blocks returns the number of divisions of the dimension length with the given // block size. func blocks(dim, bsize int) int { diff --git a/vendor/gonum.org/v1/gonum/floats/floats.go b/vendor/gonum.org/v1/gonum/floats/floats.go index 5db73a05..68c4e65c 100644 --- a/vendor/gonum.org/v1/gonum/floats/floats.go +++ b/vendor/gonum.org/v1/gonum/floats/floats.go @@ -7,6 +7,7 @@ package floats import ( "errors" "math" + "slices" "sort" "gonum.org/v1/gonum/floats/scalar" @@ -633,10 +634,10 @@ func Prod(s []float64) float64 { } // Reverse reverses the order of elements in the slice. +// +// Deprecated: This function simply calls [slices.Reverse]. func Reverse(s []float64) { - for i, j := 0, len(s)-1; i < j; i, j = i+1, j-1 { - s[i], s[j] = s[j], s[i] - } + slices.Reverse(s) } // Same returns true when the input slices have the same length and all diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dbdsqr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dbdsqr.go index e9c055b3..fd421d7e 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dbdsqr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dbdsqr.go @@ -146,7 +146,7 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v smax = math.Max(smax, math.Abs(e[i])) } - var sminl float64 + var smin float64 var thresh float64 if tol >= 0 { sminoa := math.Abs(d[0]) @@ -189,7 +189,6 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v d[m-1] = 0 } smax = math.Abs(d[m-1]) - smin := smax var l2 int var broke bool for l3 := 0; l3 < m-1; l3++ { @@ -203,7 +202,6 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v broke = true break } - smin = math.Min(smin, abss) smax = math.Max(math.Max(smax, abss), abse) } if broke { @@ -257,14 +255,14 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v if tol >= 0 { // If relative accuracy desired, apply convergence criterion forward. mu := math.Abs(d[l2]) - sminl = mu + smin = mu for l3 := l2; l3 < m-1; l3++ { if math.Abs(e[l3]) <= tol*mu { e[l3] = 0 continue Outer } mu = math.Abs(d[l3+1]) * (mu / (mu + math.Abs(e[l3]))) - sminl = math.Min(sminl, mu) + smin = math.Min(smin, mu) } } } else { @@ -277,14 +275,14 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v if tol >= 0 { // If relative accuracy desired, apply convergence criterion backward. mu := math.Abs(d[m-1]) - sminl = mu + smin = mu for l3 := m - 2; l3 >= l2; l3-- { if math.Abs(e[l3]) <= tol*mu { e[l3] = 0 continue Outer } mu = math.Abs(d[l3]) * (mu / (mu + math.Abs(e[l3]))) - sminl = math.Min(sminl, mu) + smin = math.Min(smin, mu) } } } @@ -293,7 +291,7 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v // Compute shift. First, test if shifting would ruin relative accuracy, // and if so set the shift to zero. var shift float64 - if tol >= 0 && float64(n)*tol*(sminl/smax) <= math.Max(eps, (1.0/100)*tol) { + if tol >= 0 && float64(n)*tol*(smin/smax) <= math.Max(eps, (1.0/100)*tol) { shift = 0 } else { var sl2 float64 diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgecon.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgecon.go index 1d1ca586..1d046441 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgecon.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgecon.go @@ -12,17 +12,23 @@ import ( "gonum.org/v1/gonum/lapack" ) -// Dgecon estimates the reciprocal of the condition number of the n×n matrix A -// given the LU decomposition of the matrix. The condition number computed may -// be based on the 1-norm or the ∞-norm. +// Dgecon estimates and returns the reciprocal of the condition number of the +// n×n matrix A, in either the 1-norm or the ∞-norm, using the LU factorization +// computed by Dgetrf. // -// The slice a contains the result of the LU decomposition of A as computed by Dgetrf. +// An estimate is obtained for norm(A⁻¹), and the reciprocal of the condition +// number rcond is computed as // -// anorm is the corresponding 1-norm or ∞-norm of the original matrix A. +// rcond 1 / ( norm(A) * norm(A⁻¹) ). // -// work is a temporary data slice of length at least 4*n and Dgecon will panic otherwise. +// If n is zero, rcond is always 1. // -// iwork is a temporary data slice of length at least n and Dgecon will panic otherwise. +// anorm is the 1-norm or the ∞-norm of the original matrix A. anorm must be +// non-negative, otherwise Dgecon will panic. If anorm is 0 or infinity, Dgecon +// returns 0. If anorm is NaN, Dgecon returns NaN. +// +// work must have length at least 4*n and iwork must have length at least n, +// otherwise Dgecon will panic. func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 { switch { case norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum: @@ -31,6 +37,8 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld panic(nLT0) case lda < max(1, n): panic(badLdA) + case anorm < 0: + panic(negANorm) } // Quick return if possible. @@ -48,7 +56,13 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld } // Quick return if possible. - if anorm == 0 { + switch { + case anorm == 0: + return 0 + case math.IsNaN(anorm): + // Propagate NaN. + return anorm + case math.IsInf(anorm, 1): return 0 } diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgels.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgels.go index 496e8a7e..3018973a 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgels.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgels.go @@ -131,7 +131,7 @@ func (impl Implementation) Dgels(trans blas.Transpose, m, n, nrhs int, a []float // Solve the minimization problem using a QR or an LQ decomposition. var scllen int if m >= n { - impl.Dgeqrf(m, n, a, lda, work, work[mn:], lwork-mn) + impl.Dgeqrf(m, n, a, lda, work[:n], work[mn:], lwork-mn) if trans == blas.NoTrans { impl.Dormqr(blas.Left, blas.Trans, m, nrhs, n, a, lda, diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqp3.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqp3.go index f96f03be..da8cd4fa 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqp3.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqp3.go @@ -9,35 +9,40 @@ import ( "gonum.org/v1/gonum/blas/blas64" ) -// Dgeqp3 computes a QR factorization with column pivoting of the -// m×n matrix A: A*P = Q*R using Level 3 BLAS. +// Dgeqp3 computes a QR factorization with column pivoting of the m×n matrix A: // -// The matrix Q is represented as a product of elementary reflectors +// A*P = Q*R // -// Q = H_0 H_1 . . . H_{k-1}, where k = min(m,n). +// where P is a permutation matrix, Q is an orthogonal matrix and R is a +// min(m,n)×n upper trapezoidal matrix. +// +// On return, the upper triangle of A contains the matrix R. The elements below +// the diagonal together with tau represent the matrix Q as a product of +// elementary reflectors +// +// Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n). // // Each H_i has the form // // H_i = I - tau * v * vᵀ // -// where tau and v are real vectors with v[0:i-1] = 0 and v[i] = 1; -// v[i:m] is stored on exit in A[i:m, i], and tau in tau[i]. +// where tau is a scalar and v is a vector with v[0:i] = 0 and v[i] = 1; +// v[i+1:m] is stored on exit in A[i+1:m,i], and tau in tau[i]. // -// jpvt specifies a column pivot to be applied to A. If -// jpvt[j] is at least zero, the jth column of A is permuted -// to the front of A*P (a leading column), if jpvt[j] is -1 -// the jth column of A is a free column. If jpvt[j] < -1, Dgeqp3 -// will panic. On return, jpvt holds the permutation that was -// applied; the jth column of A*P was the jpvt[j] column of A. -// jpvt must have length n or Dgeqp3 will panic. +// jpvt specifies a column pivot to be applied to A. On entry, if jpvt[j] is at +// least zero, the jth column of A is permuted to the front of A*P (a leading +// column), if jpvt[j] is -1 the jth column of A is a free column. If jpvt[j] < +// -1, Dgeqp3 will panic. On return, jpvt holds the permutation that was +// applied; the jth column of A*P was the jpvt[j] column of A. jpvt must have +// length n or Dgeqp3 will panic. // -// tau holds the scalar factors of the elementary reflectors. -// It must have length min(m, n), otherwise Dgeqp3 will panic. +// tau holds the scalar factors of the elementary reflectors. It must have +// length min(m,n), otherwise Dgeqp3 will panic. // // work must have length at least max(1,lwork), and lwork must be at least -// 3*n+1, otherwise Dgeqp3 will panic. For optimal performance lwork must -// be at least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return, -// work[0] will contain the optimal value of lwork. +// 3*n+1, otherwise Dgeqp3 will panic. For optimal performance lwork must be at +// least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return, work[0] +// will contain the optimal value of lwork. // // If lwork == -1, instead of performing Dgeqp3, only the optimal value of lwork // will be stored in work[0]. @@ -118,7 +123,7 @@ func (impl Implementation) Dgeqp3(m, n int, a []float64, lda int, jpvt []int, ta // Compute the QR factorization of nfxd columns and update remaining columns. if nfxd > 0 { na := min(m, nfxd) - impl.Dgeqrf(m, na, a, lda, tau, work, lwork) + impl.Dgeqrf(m, na, a, lda, tau[:na], work, lwork) iws = max(iws, int(work[0])) if na < n { impl.Dormqr(blas.Left, blas.Trans, m, n-na, na, a, lda, tau[:na], a[na:], lda, diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqr2.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqr2.go index c02f8a12..4d1a4b3b 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqr2.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqr2.go @@ -14,7 +14,7 @@ import "gonum.org/v1/gonum/blas" // A is modified to contain the information to construct Q and R. // The upper triangle of a contains the matrix R. The lower triangular elements // (not including the diagonal) contain the elementary reflectors. tau is modified -// to contain the reflector scales. tau must have length at least min(m,n), and +// to contain the reflector scales. tau must have length min(m,n), and // this function will panic otherwise. // // The ith elementary reflector can be explicitly constructed by first extracting @@ -57,8 +57,8 @@ func (impl Implementation) Dgeqr2(m, n int, a []float64, lda int, tau, work []fl switch { case len(a) < (m-1)*lda+n: panic(shortA) - case len(tau) < k: - panic(shortTau) + case len(tau) != k: + panic(badLenTau) } for i := 0; i < k; i++ { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqrf.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqrf.go index d14088a1..2bcbde58 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqrf.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgeqrf.go @@ -20,7 +20,7 @@ import ( // by the temporary space available. If lwork == -1, instead of performing Dgeqrf, // the optimal work length will be stored into work[0]. // -// tau must have length at least min(m,n), and this function will panic otherwise. +// tau must have length min(m,n), and this function will panic otherwise. func (impl Implementation) Dgeqrf(m, n int, a []float64, lda int, tau, work []float64, lwork int) { switch { case m < 0: @@ -52,8 +52,8 @@ func (impl Implementation) Dgeqrf(m, n int, a []float64, lda int, tau, work []fl if len(a) < (m-1)*lda+n { panic(shortA) } - if len(tau) < k { - panic(shortTau) + if len(tau) != k { + panic(badLenTau) } nbmin := 2 // Minimal block size. @@ -83,7 +83,7 @@ func (impl Implementation) Dgeqrf(m, n int, a []float64, lda int, tau, work []fl for i = 0; i < k-nx; i += nb { ib := min(k-i, nb) // Compute the QR factorization of the current block. - impl.Dgeqr2(m-i, ib, a[i*lda+i:], lda, tau[i:], work) + impl.Dgeqr2(m-i, ib, a[i*lda+i:], lda, tau[i:i+ib], work) if i+ib < n { // Form the triangular factor of the block reflector and apply Hᵀ // In Dlarft, work becomes the T matrix. diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgesvd.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgesvd.go index e0f8040f..97da749b 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgesvd.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgesvd.go @@ -406,7 +406,7 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q * R. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], work[iwork:], lwork-iwork) // Zero out below R. impl.Dlaset(blas.Lower, n-1, n-1, 0, 0, a[lda:], lda) @@ -455,14 +455,14 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float itau := ir + ldworkr*n iwork := itau + n // Compute A = Q * R. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, n, n, a, lda, work[itau:itau+n], work[iwork:], lwork-iwork) ie := itau itauq := ie + n itaup := itauq + n @@ -492,11 +492,11 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q*R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, n, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) ie := itau itauq := ie + n itaup := itauq + n @@ -537,13 +537,13 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q * R. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, n, n, a, lda, work[itau:itau+n], work[iwork:], lwork-iwork) ie := itau itauq := ie + n @@ -580,11 +580,11 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q * R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, n, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) // Copy R to VT, zeroing out below it. impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt) @@ -631,7 +631,7 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q*R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], work[iwork:], lwork-iwork) impl.Dlacpy(blas.Lower, m, n, a, lda, u, ldu) // Copy R to work[ir:], zeroing out below it. @@ -639,7 +639,7 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float 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) + impl.Dorgqr(m, m, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) ie := itau itauq := ie + n itaup := itauq + n @@ -672,11 +672,11 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q*R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, m, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) ie := itau itauq := ie + n itaup := itauq + n @@ -717,11 +717,11 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q * R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, m, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) // Copy R to work[iu:], zeroing out below it. impl.Dlacpy(blas.Upper, n, n, a, lda, work[iu:], ldworku) @@ -786,11 +786,11 @@ func (impl Implementation) Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float iwork := itau + n // Compute A = Q*R, copying result to U. - impl.Dgeqrf(m, n, a, lda, work[itau:], work[iwork:], lwork-iwork) + impl.Dgeqrf(m, n, a, lda, work[itau:itau+n], 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) + impl.Dorgqr(m, m, n, u, ldu, work[itau:itau+n], work[iwork:], lwork-iwork) // Copy R from A to VT, zeroing out below it. impl.Dlacpy(blas.Upper, n, n, a, lda, vt, ldvt) diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgetf2.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgetf2.go index b773f658..6a7003cf 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgetf2.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgetf2.go @@ -10,23 +10,27 @@ import ( "gonum.org/v1/gonum/blas/blas64" ) -// Dgetf2 computes the LU decomposition of the m×n matrix A. -// The LU decomposition is a factorization of a into +// Dgetf2 computes the LU decomposition of an m×n matrix A using partial +// pivoting with row interchanges. +// +// The LU decomposition is a factorization of A into // // A = P * L * U // -// where P is a permutation matrix, L is a unit lower triangular matrix, and -// U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored -// in place into a. +// where P is a permutation matrix, L is a lower triangular with unit diagonal +// elements (lower trapezoidal if m > n), and U is upper triangular (upper +// trapezoidal if m < n). +// +// On entry, a contains the matrix A. On return, L and U are stored in place +// into a, and P is represented by ipiv. // -// ipiv is a permutation vector. It indicates that row i of the matrix was -// changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic -// otherwise. ipiv is zero-indexed. +// ipiv contains a sequence of row interchanges. It indicates that row i of the +// matrix was interchanged with ipiv[i]. ipiv must have length min(m,n), and +// Dgetf2 will panic otherwise. ipiv is zero-indexed. // -// Dgetf2 returns whether the matrix A is singular. The LU decomposition will -// be computed regardless of the singularity of A, but division by zero -// will occur if the false is returned and the result is used to solve a -// system of equations. +// Dgetf2 returns whether the matrix A is nonsingular. The LU decomposition will +// be computed regardless of the singularity of A, but the result should not be +// used to solve a system of equation. // // Dgetf2 is an internal routine. It is exported for testing purposes. func (Implementation) Dgetf2(m, n int, a []float64, lda int, ipiv []int) (ok bool) { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgetrf.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgetrf.go index 35dde05c..38ae8efa 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dgetrf.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgetrf.go @@ -9,25 +9,27 @@ import ( "gonum.org/v1/gonum/blas/blas64" ) -// Dgetrf computes the LU decomposition of the m×n matrix A. +// Dgetrf computes the LU decomposition of an m×n matrix A using partial +// pivoting with row interchanges. +// // The LU decomposition is a factorization of A into // // A = P * L * U // -// where P is a permutation matrix, L is a unit lower triangular matrix, and -// U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored -// in place into a. +// where P is a permutation matrix, L is a lower triangular with unit diagonal +// elements (lower trapezoidal if m > n), and U is upper triangular (upper +// trapezoidal if m < n). // -// ipiv is a permutation vector. It indicates that row i of the matrix was -// changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic -// otherwise. ipiv is zero-indexed. +// On entry, a contains the matrix A. On return, L and U are stored in place +// into a, and P is represented by ipiv. // -// Dgetrf is the blocked version of the algorithm. +// ipiv contains a sequence of row interchanges. It indicates that row i of the +// matrix was interchanged with ipiv[i]. ipiv must have length min(m,n), and +// Dgetrf will panic otherwise. ipiv is zero-indexed. // -// Dgetrf returns whether the matrix A is singular. The LU decomposition will -// be computed regardless of the singularity of A, but division by zero -// will occur if the false is returned and the result is used to solve a -// system of equations. +// Dgetrf returns whether the matrix A is nonsingular. The LU decomposition will +// be computed regardless of the singularity of A, but the result should not be +// used to solve a system of equation. func (impl Implementation) Dgetrf(m, n int, a []float64, lda int, ipiv []int) (ok bool) { mn := min(m, n) switch { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dgghrd.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dgghrd.go new file mode 100644 index 00000000..c9d6b4d1 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dgghrd.go @@ -0,0 +1,125 @@ +// Copyright ©2023 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 ( + "gonum.org/v1/gonum/blas" + "gonum.org/v1/gonum/blas/blas64" + "gonum.org/v1/gonum/lapack" +) + +// Dgghrd reduces a pair of real matrices (A,B) to generalized upper Hessenberg +// form using orthogonal transformations, where A is a general matrix and B is +// upper triangular. +// +// This subroutine simultaneously reduces A to a Hessenberg matrix H +// +// Qᵀ*A*Z = H, +// +// and transforms B to another upper triangular matrix T +// +// Qᵀ*B*Z = T. +// +// The orthogonal matrices Q and Z are determined as products of Givens +// rotations. They may either be formed explicitly (lapack.OrthoExplicit), or +// they may be postmultiplied into input matrices Q1 and Z1 +// (lapack.OrthoPostmul), so that +// +// Q1 * A * Z1ᵀ = (Q1*Q) * H * (Z1*Z)ᵀ, +// Q1 * B * Z1ᵀ = (Q1*Q) * T * (Z1*Z)ᵀ. +// +// ilo and ihi determine the block of A that will be reduced. It must hold that +// +// - 0 <= ilo <= ihi < n if n > 0, +// - ilo == 0 and ihi == -1 if n == 0, +// +// otherwise Dgghrd will panic. +// +// Dgghrd is an internal routine. It is exported for testing purposes. +func (impl Implementation) Dgghrd(compq, compz lapack.OrthoComp, n, ilo, ihi int, a []float64, lda int, b []float64, ldb int, q []float64, ldq int, z []float64, ldz int) { + switch { + case compq != lapack.OrthoNone && compq != lapack.OrthoExplicit && compq != lapack.OrthoPostmul: + panic(badOrthoComp) + case compz != lapack.OrthoNone && compz != lapack.OrthoExplicit && compz != lapack.OrthoPostmul: + panic(badOrthoComp) + case n < 0: + panic(nLT0) + case ilo < 0 || max(0, n-1) < ilo: + panic(badIlo) + case ihi < min(ilo, n-1) || n <= ihi: + panic(badIhi) + case lda < max(1, n): + panic(badLdA) + case ldb < max(1, n): + panic(badLdB) + case (compq != lapack.OrthoNone && ldq < n) || ldq < 1: + panic(badLdQ) + case (compz != lapack.OrthoNone && ldz < n) || ldz < 1: + panic(badLdZ) + } + + // Quick return if possible. + if n == 0 { + return + } + + switch { + case len(a) < (n-1)*lda+n: + panic(shortA) + case len(b) < (n-1)*ldb+n: + panic(shortB) + case compq != lapack.OrthoNone && len(q) < (n-1)*ldq+n: + panic(shortQ) + case compz != lapack.OrthoNone && len(z) < (n-1)*ldz+n: + panic(shortZ) + } + + if compq == lapack.OrthoExplicit { + impl.Dlaset(blas.All, n, n, 0, 1, q, ldq) + } + if compz == lapack.OrthoExplicit { + impl.Dlaset(blas.All, n, n, 0, 1, z, ldz) + } + + // Quick return if possible. + if n == 1 { + return + } + + // Zero out lower triangle of B. + for i := 1; i < n; i++ { + for j := 0; j < i; j++ { + b[i*ldb+j] = 0 + } + } + bi := blas64.Implementation() + // Reduce A and B. + for jcol := ilo; jcol <= ihi-2; jcol++ { + for jrow := ihi; jrow >= jcol+2; jrow-- { + // Step 1: rotate rows jrow-1, jrow to kill A[jrow,jcol]. + var c, s float64 + c, s, a[(jrow-1)*lda+jcol] = impl.Dlartg(a[(jrow-1)*lda+jcol], a[jrow*lda+jcol]) + a[jrow*lda+jcol] = 0 + + bi.Drot(n-jcol-1, a[(jrow-1)*lda+jcol+1:], 1, a[jrow*lda+jcol+1:], 1, c, s) + bi.Drot(n+2-jrow-1, b[(jrow-1)*ldb+jrow-1:], 1, b[jrow*ldb+jrow-1:], 1, c, s) + + if compq != lapack.OrthoNone { + bi.Drot(n, q[jrow-1:], ldq, q[jrow:], ldq, c, s) + } + + // Step 2: rotate columns jrow, jrow-1 to kill B[jrow,jrow-1]. + c, s, b[jrow*ldb+jrow] = impl.Dlartg(b[jrow*ldb+jrow], b[jrow*ldb+jrow-1]) + b[jrow*ldb+jrow-1] = 0 + + bi.Drot(ihi+1, a[jrow:], lda, a[jrow-1:], lda, c, s) + bi.Drot(jrow, b[jrow:], ldb, b[jrow-1:], ldb, c, s) + + if compz != lapack.OrthoNone { + bi.Drot(n, z[jrow:], ldz, z[jrow-1:], ldz, c, s) + } + } + } +} diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dggsvp3.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dggsvp3.go index ace13c6c..f7f04c76 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dggsvp3.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dggsvp3.go @@ -154,7 +154,7 @@ func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, if p > 1 { impl.Dlacpy(blas.Lower, p-1, min(p, n), b[ldb:], ldb, v[ldv:], ldv) } - impl.Dorg2r(p, p, min(p, n), v, ldv, tau, work) + impl.Dorg2r(p, p, min(p, n), v, ldv, tau[:min(p, n)], work) } // Clean up B. @@ -216,7 +216,7 @@ func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, } // Update A12 := Uᵀ*A12, where A12 = A[0:m, n-l:n]. - impl.Dorm2r(blas.Left, blas.Trans, m, l, min(m, n-l), a, lda, tau, a[n-l:], lda, work) + impl.Dorm2r(blas.Left, blas.Trans, m, l, min(m, n-l), a, lda, tau[:min(m, n-l)], a[n-l:], lda, work) if wantu { // Copy the details of U, and form U. @@ -224,7 +224,8 @@ func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, if m > 1 { impl.Dlacpy(blas.Lower, m-1, min(m, n-l), a[lda:], lda, u[ldu:], ldu) } - impl.Dorg2r(m, m, min(m, n-l), u, ldu, tau, work) + k := min(m, n-l) + impl.Dorg2r(m, m, k, u, ldu, tau[:k], work) } if wantq { @@ -250,7 +251,7 @@ func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, if wantq { // Update Q[0:n, 0:n-l] := Q[0:n, 0:n-l]*Z1ᵀ. - impl.Dorm2r(blas.Right, blas.Trans, n, n-l, k, a, lda, tau, q, ldq, work) + impl.Dorm2r(blas.Right, blas.Trans, n, n-l, k, a, lda, tau[:k], q, ldq, work) } // Clean up A. @@ -265,10 +266,10 @@ func (impl Implementation) Dggsvp3(jobU, jobV, jobQ lapack.GSVDJob, m, p, n int, if m > k { // QR factorization of A[k:m, n-l:n]. - impl.Dgeqr2(m-k, l, a[k*lda+n-l:], lda, tau, work) + impl.Dgeqr2(m-k, l, a[k*lda+n-l:], lda, tau[:min(m-k, l)], work) if wantu { // Update U[:, k:m) := U[:, k:m]*U1. - impl.Dorm2r(blas.Right, blas.NoTrans, m, m-k, min(m-k, l), a[k*lda+n-l:], lda, tau, u[k:], ldu, work) + impl.Dorm2r(blas.Right, blas.NoTrans, m, m-k, min(m-k, l), a[k*lda+n-l:], lda, tau[:min(m-k, l)], u[k:], ldu, work) } // Clean up A. diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dhseqr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dhseqr.go index 80fe19bb..beccf132 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dhseqr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dhseqr.go @@ -191,7 +191,7 @@ func (impl Implementation) Dhseqr(job lapack.SchurJob, compz lapack.SchurComp, n // Matrices of order ntiny or smaller must be processed by // Dlahqr because of insufficient subdiagonal scratch space. // This is a hard limit. - ntiny = 11 + ntiny = 15 // nl is the size of a local workspace to help small matrices // through a rare Dlahqr failure. nl > ntiny is required and diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlahqr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlahqr.go index 13f28560..6f120254 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlahqr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlahqr.go @@ -150,6 +150,9 @@ func (impl Implementation) Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float6 itmax := 30 * max(10, nh) // Total number of QR iterations allowed. + // kdefl counts the number of iterations since a deflation. + kdefl := 0 + // The main loop begins here. i is the loop index and decreases from ihi // to ilo in steps of 1 or 2. Each iteration of the loop works with the // active submatrix in rows and columns l to i. Eigenvalues i+1 to ihi @@ -207,6 +210,7 @@ func (impl Implementation) Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float6 converged = true break } + kdefl++ // Now the active submatrix is in rows and columns l to // i. If eigenvalues only are being computed, only the @@ -217,20 +221,21 @@ func (impl Implementation) Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float6 } const ( - dat1 = 3.0 - dat2 = -0.4375 + dat1 = 0.75 + dat2 = -0.4375 + kexsh = 10 ) var h11, h21, h12, h22 float64 - switch its { - case 10: // Exceptional shift. - s := math.Abs(h[(l+1)*ldh+l]) + math.Abs(h[(l+2)*ldh+l+1]) - h11 = dat1*s + h[l*ldh+l] + switch { + case kdefl%(2*kexsh) == 0: // Exceptional shift. + s := math.Abs(h[i*ldh+i-1]) + math.Abs(h[(i-1)*ldh+i-2]) + h11 = dat1*s + h[i*ldh+i] h12 = dat2 * s h21 = s h22 = h11 - case 20: // Exceptional shift. - s := math.Abs(h[i*ldh+i-1]) + math.Abs(h[(i-1)*ldh+i-2]) - h11 = dat1*s + h[i*ldh+i] + case kdefl%kexsh == 0: // Exceptional shift. + s := math.Abs(h[(l+1)*ldh+l]) + math.Abs(h[(l+2)*ldh+l+1]) + h11 = dat1*s + h[l*ldh+l] h12 = dat2 * s h21 = s h22 = h11 @@ -434,6 +439,9 @@ func (impl Implementation) Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float6 } } + // Reset deflation counter. + kdefl = 0 + // Return to start of the main loop with new value of i. i = l - 1 } diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlanhs.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlanhs.go new file mode 100644 index 00000000..054b90f0 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlanhs.go @@ -0,0 +1,78 @@ +// Copyright ©2023 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/blas64" + "gonum.org/v1/gonum/lapack" +) + +// Dlanhs returns the value of the one norm, or the Frobenius norm, or the +// infinity norm, or the element of largest absolute value of a Hessenberg +// matrix A. +// +// If norm is lapack.MaxColumnSum, work must have length at least n. +func (impl Implementation) Dlanhs(norm lapack.MatrixNorm, n int, a []float64, lda int, work []float64) float64 { + switch { + case norm != lapack.MaxRowSum && norm != lapack.MaxAbs && norm != lapack.MaxColumnSum && norm != lapack.Frobenius: + panic(badNorm) + case n < 0: + panic(nLT0) + case lda < max(1, n): + panic(badLdA) + } + + if n == 0 { + return 0 + } + + switch { + case len(a) < (n-1)*lda+n: + panic(shortA) + case norm == lapack.MaxColumnSum && len(work) < n: + panic(shortWork) + } + + bi := blas64.Implementation() + var value float64 + switch norm { + case lapack.MaxAbs: + for i := 0; i < n; i++ { + minj := max(0, i-1) + for _, v := range a[i*lda+minj : i*lda+n] { + value = math.Max(value, math.Abs(v)) + } + } + case lapack.MaxColumnSum: + for i := 0; i < n; i++ { + work[i] = 0 + } + for i := 0; i < n; i++ { + for j := max(0, i-1); j < n; j++ { + work[j] += math.Abs(a[i*lda+j]) + } + } + for _, v := range work[:n] { + value = math.Max(value, v) + } + case lapack.MaxRowSum: + for i := 0; i < n; i++ { + minj := max(0, i-1) + sum := bi.Dasum(n-minj, a[i*lda+minj:], 1) + value = math.Max(value, sum) + } + case lapack.Frobenius: + scale := 0.0 + sum := 1.0 + for i := 0; i < n; i++ { + minj := max(0, i-1) + scale, sum = impl.Dlassq(n-minj, a[i*lda+minj:], 1, scale, sum) + } + value = scale * math.Sqrt(sum) + } + return value +} diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr04.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr04.go index 3faaa2fc..8e4b266b 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr04.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr04.go @@ -118,7 +118,7 @@ func (impl Implementation) Dlaqr04(wantt, wantz bool, n, ilo, ihi int, h []float // Matrices of order ntiny or smaller must be processed by // Dlahqr because of insufficient subdiagonal scratch space. // This is a hard limit. - ntiny = 11 + ntiny = 15 // Exceptional deflation windows: try to cure rare slow // convergence by varying the size of the deflation window after // kexnw iterations. @@ -218,20 +218,20 @@ func (impl Implementation) Dlaqr04(wantt, wantz bool, n, ilo, ihi int, h []float } else { fname = "DLAQR4" } - // nwr is the recommended deflation window size. n is greater than 11, + // nwr is the recommended deflation window size. n is greater than ntiny, // so there is enough subdiagonal workspace for nwr >= 2 as required. - // (In fact, there is enough subdiagonal space for nwr >= 3.) - // TODO(vladimir-ch): If there is enough space for nwr >= 3, should we + // (In fact, there is enough subdiagonal space for nwr >= 4.) + // TODO(vladimir-ch): If there is enough space for nwr >= 4, should we // use it? nwr := impl.Ilaenv(13, fname, jbcmpz, n, ilo, ihi, lwork) nwr = max(2, nwr) nwr = min(ihi-ilo+1, min((n-1)/3, nwr)) - // nsr is the recommended number of simultaneous shifts. n is greater - // than 11, so there is enough subdiagonal workspace for nsr to be even - // and greater than or equal to two as required. + // nsr is the recommended number of simultaneous shifts. n is greater than + // ntiny, so there is enough subdiagonal workspace for nsr to be even and + // greater than or equal to two as required. nsr := impl.Ilaenv(15, fname, jbcmpz, n, ilo, ihi, lwork) - nsr = min(nsr, min((n+6)/9, ihi-ilo)) + nsr = min(nsr, min((n-3)/6, ihi-ilo)) nsr = max(2, nsr&^1) // Workspace query call to Dlaqr23. @@ -264,7 +264,7 @@ func (impl Implementation) Dlaqr04(wantt, wantz bool, n, ilo, ihi int, h []float // nsmax is the largest number of simultaneous shifts for which there is // sufficient workspace. - nsmax := min((n+6)/9, 2*lwork/3) &^ 1 + nsmax := min((n-3)/6, 2*lwork/3) &^ 1 ndfl := 1 // Number of iterations since last deflation. ndec := 0 // Deflation window size decrement. @@ -470,7 +470,7 @@ func (impl Implementation) Dlaqr04(wantt, wantz bool, n, ilo, ihi int, h []float // (nhv must be at least kdu but more is better), // - an nhv×kdu vertical work array WV along the left-hand-edge // (nhv must be at least kdu but more is better). - kdu := 3*ns - 3 + kdu := 2 * ns ku := n - kdu kwh := kdu kwv = kdu + 3 diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr5.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr5.go index 43b425b8..443a53d5 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr5.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr5.go @@ -24,9 +24,8 @@ import ( // multiply to update far-from-diagonal matrix entries. // 1: Dlaqr5 will accumulate reflections and use matrix-matrix multiply to // update far-from-diagonal matrix entries. -// 2: Dlaqr5 will accumulate reflections, use matrix-matrix multiply to update -// far-from-diagonal matrix entries, and take advantage of 2×2 block -// structure during matrix multiplies. +// 2: Same as kacc22=1. This option used to enable exploiting the 2×2 structure +// during matrix multiplications, but this is no longer supported. // // For other values of kacc2 Dlaqr5 will panic. // @@ -63,11 +62,11 @@ import ( // v and ldv represent an auxiliary matrix V of size (nshfts/2)×3. Note that V // is transposed with respect to the reference netlib implementation. // -// u and ldu represent an auxiliary matrix of size (3*nshfts-3)×(3*nshfts-3). +// u and ldu represent an auxiliary matrix of size (2*nshfts)×(2*nshfts). // -// wh and ldwh represent an auxiliary matrix of size (3*nshfts-3)×nh. +// wh and ldwh represent an auxiliary matrix of size (2*nshfts-1)×nh. // -// wv and ldwv represent an auxiliary matrix of size nv×(3*nshfts-3). +// wv and ldwv represent an auxiliary matrix of size nv×(2*nshfts-1). // // Dlaqr5 is an internal routine. It is exported for testing purposes. func (impl Implementation) Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, nshfts int, sr, si []float64, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, v []float64, ldv int, u []float64, ldu int, nv int, wv []float64, ldwv int, nh int, wh []float64, ldwh int) { @@ -110,23 +109,23 @@ func (impl Implementation) Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, case len(v) < (nshfts/2-1)*ldv+3: panic(shortV) - case ldu < max(1, 3*nshfts-3): + case ldu < max(1, 2*nshfts): panic(badLdU) - case len(u) < (3*nshfts-3-1)*ldu+3*nshfts-3: + case len(u) < (2*nshfts-1)*ldu+2*nshfts: panic(shortU) case nv < 0: panic(nvLT0) - case ldwv < max(1, 3*nshfts-3): + case ldwv < max(1, 2*nshfts): panic(badLdWV) - case len(wv) < (nv-1)*ldwv+3*nshfts-3: + case len(wv) < (nv-1)*ldwv+2*nshfts: panic(shortWV) case nh < 0: panic(nhLT0) case ldwh < max(1, nh): panic(badLdWH) - case len(wh) < (3*nshfts-3-1)*ldwh+nh: + case len(wh) < (2*nshfts-1)*ldwh+nh: panic(shortWH) case ktop > 0 && h[ktop*ldh+ktop-1] != 0: @@ -174,8 +173,6 @@ func (impl Implementation) Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, // Use accumulated reflections to update far-from-diagonal entries? accum := kacc22 == 1 || kacc22 == 2 - // If so, exploit the 2×2 block structure? - blk22 := ns > 2 && kacc22 == 2 // Clear trash. if ktop+2 <= kbot { @@ -186,103 +183,49 @@ func (impl Implementation) Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, nbmps := ns / 2 // kdu = width of slab. - kdu := 6*nbmps - 3 + kdu := 4 * nbmps // Create and chase chains of nbmps bulges. - for incol := 3*(1-nbmps) + ktop - 1; incol <= kbot-2; incol += 3*nbmps - 2 { + for incol := ktop - 2*nbmps + 1; incol <= kbot-2; incol += 2 * nbmps { + // jtop is an index from which updates from the right start. + var jtop int + switch { + case accum: + jtop = max(ktop, incol) + case wantt: + default: + jtop = ktop + } ndcol := incol + kdu if accum { impl.Dlaset(blas.All, kdu, kdu, 0, 1, u, ldu) } - // Near-the-diagonal bulge chase. The following loop performs // the near-the-diagonal part of a small bulge multi-shift QR - // sweep. Each 6*nbmps-2 column diagonal chunk extends from + // sweep. Each 4*nbmps column diagonal chunk extends from // column incol to column ndcol (including both column incol and - // column ndcol). The following loop chases a 3*nbmps column - // long chain of nbmps bulges 3*nbmps-2 columns to the right. + // column ndcol). The following loop chases a 2*nbmps+1 column + // long chain of nbmps bulges 2*nbmps columns to the right. // (incol may be less than ktop and ndcol may be greater than // kbot indicating phantom columns from which to chase bulges // before they are actually introduced or to which to chase // bulges beyond column kbot.) - for krcol := incol; krcol <= min(incol+3*nbmps-3, kbot-2); krcol++ { + for krcol := incol; krcol <= min(incol+2*nbmps-1, kbot-2); krcol++ { // Bulges number mtop to mbot are active double implicit // shift bulges. There may or may not also be small 2×2 // bulge, if there is room. The inactive bulges (if any) // must wait until the active bulges have moved down the // diagonal to make room. The phantom matrix paradigm // described above helps keep track. - - mtop := max(0, ((ktop-1)-krcol+2)/3) - mbot := min(nbmps, (kbot-krcol)/3) - 1 + mtop := max(0, (ktop-krcol)/2) + mbot := min(nbmps, (kbot-krcol-1)/2) - 1 m22 := mbot + 1 - bmp22 := (mbot < nbmps-1) && (krcol+3*m22 == kbot-2) - - // Generate reflections to chase the chain right one - // column. (The minimum value of k is ktop-1.) - for m := mtop; m <= mbot; m++ { - k := krcol + 3*m - if k == ktop-1 { - impl.Dlaqr1(3, h[ktop*ldh+ktop:], ldh, - sr[2*m], si[2*m], sr[2*m+1], si[2*m+1], - v[m*ldv:m*ldv+3]) - alpha := v[m*ldv] - _, v[m*ldv] = impl.Dlarfg(3, alpha, v[m*ldv+1:m*ldv+3], 1) - continue - } - beta := h[(k+1)*ldh+k] - v[m*ldv+1] = h[(k+2)*ldh+k] - v[m*ldv+2] = h[(k+3)*ldh+k] - beta, v[m*ldv] = impl.Dlarfg(3, beta, v[m*ldv+1:m*ldv+3], 1) - - // A bulge may collapse because of vigilant deflation or - // destructive underflow. In the underflow case, try the - // two-small-subdiagonals trick to try to reinflate the - // bulge. - if h[(k+3)*ldh+k] != 0 || h[(k+3)*ldh+k+1] != 0 || h[(k+3)*ldh+k+2] == 0 { - // Typical case: not collapsed (yet). - h[(k+1)*ldh+k] = beta - h[(k+2)*ldh+k] = 0 - h[(k+3)*ldh+k] = 0 - continue - } - - // Atypical case: collapsed. Attempt to reintroduce - // ignoring H[k+1,k] and H[k+2,k]. If the fill - // resulting from the new reflector is too large, - // then abandon it. Otherwise, use the new one. - var vt [3]float64 - impl.Dlaqr1(3, h[(k+1)*ldh+k+1:], ldh, sr[2*m], - si[2*m], sr[2*m+1], si[2*m+1], vt[:]) - alpha := vt[0] - _, vt[0] = impl.Dlarfg(3, alpha, vt[1:3], 1) - refsum := vt[0] * (h[(k+1)*ldh+k] + vt[1]*h[(k+2)*ldh+k]) - - dsum := math.Abs(h[k*ldh+k]) + math.Abs(h[(k+1)*ldh+k+1]) + math.Abs(h[(k+2)*ldh+k+2]) - if math.Abs(h[(k+2)*ldh+k]-refsum*vt[1])+math.Abs(refsum*vt[2]) > ulp*dsum { - // Starting a new bulge here would create - // non-negligible fill. Use the old one with - // trepidation. - h[(k+1)*ldh+k] = beta - h[(k+2)*ldh+k] = 0 - h[(k+3)*ldh+k] = 0 - continue - } else { - // Starting a new bulge here would create - // only negligible fill. Replace the old - // reflector with the new one. - h[(k+1)*ldh+k] -= refsum - h[(k+2)*ldh+k] = 0 - h[(k+3)*ldh+k] = 0 - v[m*ldv] = vt[0] - v[m*ldv+1] = vt[1] - v[m*ldv+2] = vt[2] - } - } - - // Generate a 2×2 reflection, if needed. + bmp22 := (mbot < nbmps-1) && (krcol+2*m22 == kbot-2) + // Generate reflections to chase the chain right one column. + // The minimum value of k is ktop-1. if bmp22 { - k := krcol + 3*m22 + // Special case: 2×2 reflection at bottom treated separately. + k := krcol + 2*m22 if k == ktop-1 { impl.Dlaqr1(2, h[(k+1)*ldh+k+1:], ldh, sr[2*m22], si[2*m22], sr[2*m22+1], si[2*m22+1], @@ -296,353 +239,322 @@ func (impl Implementation) Dlaqr5(wantt, wantz bool, kacc22 int, n, ktop, kbot, h[(k+1)*ldh+k] = beta h[(k+2)*ldh+k] = 0 } - } - - // Multiply H by reflections from the left. - var jbot int - switch { - case accum: - jbot = min(ndcol, kbot) - case wantt: - jbot = n - 1 - default: - jbot = kbot - } - for j := max(ktop, krcol); j <= jbot; j++ { - mend := min(mbot+1, (j-krcol+2)/3) - 1 - for m := mtop; m <= mend; m++ { - k := krcol + 3*m - refsum := v[m*ldv] * (h[(k+1)*ldh+j] + - v[m*ldv+1]*h[(k+2)*ldh+j] + v[m*ldv+2]*h[(k+3)*ldh+j]) - h[(k+1)*ldh+j] -= refsum - h[(k+2)*ldh+j] -= refsum * v[m*ldv+1] - h[(k+3)*ldh+j] -= refsum * v[m*ldv+2] - } - } - if bmp22 { - k := krcol + 3*m22 - for j := max(k+1, ktop); j <= jbot; j++ { - refsum := v[m22*ldv] * (h[(k+1)*ldh+j] + v[m22*ldv+1]*h[(k+2)*ldh+j]) - h[(k+1)*ldh+j] -= refsum - h[(k+2)*ldh+j] -= refsum * v[m22*ldv+1] + // Perform update from right within computational window. + t1 := v[m22*ldv] + t2 := t1 * v[m22*ldv+1] + for j := jtop; j <= min(kbot, k+3); j++ { + refsum := h[j*ldh+k+1] + v[m22*ldv+1]*h[j*ldh+k+2] + h[j*ldh+k+1] -= refsum * t1 + h[j*ldh+k+2] -= refsum * t2 } - } - - // Multiply H by reflections from the right. Delay filling in the last row - // until the vigilant deflation check is complete. - var jtop int - switch { - case accum: - jtop = max(ktop, incol) - case wantt: - jtop = 0 - default: - jtop = ktop - } - for m := mtop; m <= mbot; m++ { - if v[m*ldv] == 0 { - continue + // Perform update from left within computational window. + var jbot int + switch { + case accum: + jbot = min(ndcol, kbot) + case wantt: + jbot = n - 1 + default: + jbot = kbot } - k := krcol + 3*m - for j := jtop; j <= min(kbot, k+3); j++ { - refsum := v[m*ldv] * (h[j*ldh+k+1] + - v[m*ldv+1]*h[j*ldh+k+2] + v[m*ldv+2]*h[j*ldh+k+3]) - h[j*ldh+k+1] -= refsum - h[j*ldh+k+2] -= refsum * v[m*ldv+1] - h[j*ldh+k+3] -= refsum * v[m*ldv+2] + t1 = v[m22*ldv] + t2 = t1 * v[m22*ldv+1] + for j := k + 1; j <= jbot; j++ { + refsum := h[(k+1)*ldh+j] + v[m22*ldv+1]*h[(k+2)*ldh+j] + h[(k+1)*ldh+j] -= refsum * t1 + h[(k+2)*ldh+j] -= refsum * t2 } - if accum { - // Accumulate U. (If necessary, update Z later with an - // efficient matrix-matrix multiply.) - kms := k - incol - for j := max(0, ktop-incol-1); j < kdu; j++ { - refsum := v[m*ldv] * (u[j*ldu+kms] + - v[m*ldv+1]*u[j*ldu+kms+1] + v[m*ldv+2]*u[j*ldu+kms+2]) - u[j*ldu+kms] -= refsum - u[j*ldu+kms+1] -= refsum * v[m*ldv+1] - u[j*ldu+kms+2] -= refsum * v[m*ldv+2] + // The following convergence test requires that the traditional + // small-compared-to-nearby-diagonals criterion and the Ahues & + // Tisseur (LAWN 122, 1997) criteria both be satisfied. The latter + // improves accuracy in some examples. Falling back on an alternate + // convergence criterion when tst1 or tst2 is zero (as done here) is + // traditional but probably unnecessary. + if k >= ktop && h[(k+1)*ldh+k] != 0 { + tst1 := math.Abs(h[k*ldh+k]) + math.Abs(h[(k+1)*ldh+k+1]) + if tst1 == 0 { + if k >= ktop+1 { + tst1 += math.Abs(h[k*ldh+k-1]) + } + if k >= ktop+2 { + tst1 += math.Abs(h[k*ldh+k-2]) + } + if k >= ktop+3 { + tst1 += math.Abs(h[k*ldh+k-3]) + } + if k <= kbot-2 { + tst1 += math.Abs(h[(k+2)*ldh+k+1]) + } + if k <= kbot-3 { + tst1 += math.Abs(h[(k+3)*ldh+k+1]) + } + if k <= kbot-4 { + tst1 += math.Abs(h[(k+4)*ldh+k+1]) + } } - } else if wantz { - // U is not accumulated, so update Z now by multiplying by - // reflections from the right. - for j := iloz; j <= ihiz; j++ { - refsum := v[m*ldv] * (z[j*ldz+k+1] + - v[m*ldv+1]*z[j*ldz+k+2] + v[m*ldv+2]*z[j*ldz+k+3]) - z[j*ldz+k+1] -= refsum - z[j*ldz+k+2] -= refsum * v[m*ldv+1] - z[j*ldz+k+3] -= refsum * v[m*ldv+2] + if math.Abs(h[(k+1)*ldh+k]) <= math.Max(smlnum, ulp*tst1) { + h12 := math.Max(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) + h21 := math.Min(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) + h11 := math.Max(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) + h22 := math.Min(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) + scl := h11 + h12 + tst2 := h22 * (h11 / scl) + if tst2 == 0 || h21*(h12/scl) <= math.Max(smlnum, ulp*tst2) { + h[(k+1)*ldh+k] = 0 + } } } - } - - // Special case: 2×2 reflection (if needed). - if bmp22 && v[m22*ldv] != 0 { - k := krcol + 3*m22 - for j := jtop; j <= min(kbot, k+3); j++ { - refsum := v[m22*ldv] * (h[j*ldh+k+1] + v[m22*ldv+1]*h[j*ldh+k+2]) - h[j*ldh+k+1] -= refsum - h[j*ldh+k+2] -= refsum * v[m22*ldv+1] - } + // Accumulate orthogonal transformations. if accum { - kms := k - incol + kms := k - incol - 1 + t1 = v[m22*ldv] + t2 = t1 * v[m22*ldv+1] for j := max(0, ktop-incol-1); j < kdu; j++ { - refsum := v[m22*ldv] * (u[j*ldu+kms] + v[m22*ldv+1]*u[j*ldu+kms+1]) - u[j*ldu+kms] -= refsum - u[j*ldu+kms+1] -= refsum * v[m22*ldv+1] + refsum := u[j*ldu+kms+1] + v[m22*ldv+1]*u[j*ldu+kms+2] + u[j*ldu+kms+1] -= refsum * t1 + u[j*ldu+kms+2] -= refsum * t2 } } else if wantz { + t1 = v[m22*ldv] + t2 = t1 * v[m22*ldv+1] for j := iloz; j <= ihiz; j++ { - refsum := v[m22*ldv] * (z[j*ldz+k+1] + v[m22*ldv+1]*z[j*ldz+k+2]) - z[j*ldz+k+1] -= refsum - z[j*ldz+k+2] -= refsum * v[m22*ldv+1] + refsum := z[j*ldz+k+1] + v[m22*ldv+1]*z[j*ldz+k+2] + z[j*ldz+k+1] -= refsum * t1 + z[j*ldz+k+2] -= refsum * t2 } } } - - // Vigilant deflation check. - mstart := mtop - if krcol+3*mstart < ktop { - mstart++ - } - mend := mbot - if bmp22 { - mend++ - } - if krcol == kbot-2 { - mend++ - } - for m := mstart; m <= mend; m++ { - k := min(kbot-1, krcol+3*m) - + // Normal case: Chain of 3×3 reflections. + for m := mbot; m >= mtop; m-- { + k := krcol + 2*m + if k == ktop-1 { + impl.Dlaqr1(3, h[ktop*ldh+ktop:], ldh, + sr[2*m], si[2*m], sr[2*m+1], si[2*m+1], + v[m*ldv:m*ldv+3]) + alpha := v[m*ldv] + _, v[m*ldv] = impl.Dlarfg(3, alpha, v[m*ldv+1:m*ldv+3], 1) + } else { + // Perform delayed transformation of row below m-th bulge. + // Exploit fact that first two elements of row are actually + // zero. + t1 := v[m*ldv] + t2 := t1 * v[m*ldv+1] + t3 := t1 * v[m*ldv+2] + refsum := v[m*ldv+2] * h[(k+3)*ldh+k+2] + h[(k+3)*ldh+k] = -refsum * t1 + h[(k+3)*ldh+k+1] = -refsum * t2 + h[(k+3)*ldh+k+2] -= refsum * t3 + // Calculate reflection to move m-th bulge one step. + beta := h[(k+1)*ldh+k] + v[m*ldv+1] = h[(k+2)*ldh+k] + v[m*ldv+2] = h[(k+3)*ldh+k] + beta, v[m*ldv] = impl.Dlarfg(3, beta, v[m*ldv+1:m*ldv+3], 1) + // A bulge may collapse because of vigilant deflation or + // destructive underflow. In the underflow case, try the + // two-small-subdiagonals trick to try to reinflate the + // bulge. + if h[(k+3)*ldh+k] != 0 || h[(k+3)*ldh+k+1] != 0 || h[(k+3)*ldh+k+2] == 0 { + // Typical case: not collapsed (yet). + h[(k+1)*ldh+k] = beta + h[(k+2)*ldh+k] = 0 + h[(k+3)*ldh+k] = 0 + } else { + // Atypical case: collapsed. Attempt to reintroduce + // ignoring H[k+1,k] and H[k+2,k]. If the fill resulting + // from the new reflector is too large, then abandon it. + // Otherwise, use the new one. + var vt [3]float64 + impl.Dlaqr1(3, h[(k+1)*ldh+k+1:], ldh, + sr[2*m], si[2*m], sr[2*m+1], si[2*m+1], + vt[:]) + _, vt[0] = impl.Dlarfg(3, vt[0], vt[1:3], 1) + t1 = vt[0] + t2 = t1 * vt[1] + t3 = t1 * vt[2] + refsum = h[(k+1)*ldh+k] + vt[1]*h[(k+2)*ldh+k] + dsum := math.Abs(h[k*ldh+k]) + math.Abs(h[(k+1)*ldh+k+1]) + math.Abs(h[(k+2)*ldh+k+2]) + if math.Abs(h[(k+2)*ldh+k]-refsum*t2)+math.Abs(refsum*t3) > ulp*dsum { + // Starting a new bulge here would create + // non-negligible fill. Use the old one with + // trepidation. + h[(k+1)*ldh+k] = beta + h[(k+2)*ldh+k] = 0 + h[(k+3)*ldh+k] = 0 + } else { + // Starting a new bulge here would create only + // negligible fill. Replace the old reflector with + // the new one. + h[(k+1)*ldh+k] -= refsum * t1 + h[(k+2)*ldh+k] = 0 + h[(k+3)*ldh+k] = 0 + v[m*ldv] = vt[0] + v[m*ldv+1] = vt[1] + v[m*ldv+2] = vt[2] + } + } + } + // Apply reflection from the right and the first column of + // update from the left. These updates are required for the + // vigilant deflation check. We still delay most of the updates + // from the left for efficiency. + t1 := v[m*ldv] + t2 := t1 * v[m*ldv+1] + t3 := t1 * v[m*ldv+2] + for j := jtop; j <= min(kbot, k+3); j++ { + refsum := h[j*ldh+k+1] + v[m*ldv+1]*h[j*ldh+k+2] + v[m*ldv+2]*h[j*ldh+k+3] + h[j*ldh+k+1] -= refsum * t1 + h[j*ldh+k+2] -= refsum * t2 + h[j*ldh+k+3] -= refsum * t3 + } + // Perform update from left for subsequent column. + refsum := h[(k+1)*ldh+k+1] + v[m*ldv+1]*h[(k+2)*ldh+k+1] + v[m*ldv+2]*h[(k+3)*ldh+k+1] + h[(k+1)*ldh+k+1] -= refsum * t1 + h[(k+2)*ldh+k+1] -= refsum * t2 + h[(k+3)*ldh+k+1] -= refsum * t3 // The following convergence test requires that the tradition // small-compared-to-nearby-diagonals criterion and the Ahues & - // Tisseur (LAWN 122, 1997) criteria both be satisfied. The latter - // improves accuracy in some examples. Falling back on an alternate - // convergence criterion when tst1 or tst2 is zero (as done here) is - // traditional but probably unnecessary. - - if h[(k+1)*ldh+k] == 0 { + // Tisseur (LAWN 122, 1997) criteria both be satisfied. The + // latter improves accuracy in some examples. Falling back on an + // alternate convergence criterion when tst1 or tst2 is zero (as + // done here) is traditional but probably unnecessary. + if k < ktop { continue } - tst1 := math.Abs(h[k*ldh+k]) + math.Abs(h[(k+1)*ldh+k+1]) - if tst1 == 0 { - if k >= ktop+1 { - tst1 += math.Abs(h[k*ldh+k-1]) - } - if k >= ktop+2 { - tst1 += math.Abs(h[k*ldh+k-2]) - } - if k >= ktop+3 { - tst1 += math.Abs(h[k*ldh+k-3]) + if h[(k+1)*ldh+k] != 0 { + tst1 := math.Abs(h[k*ldh+k]) + math.Abs(h[(k+1)*ldh+k+1]) + if tst1 == 0 { + if k >= ktop+1 { + tst1 += math.Abs(h[k*ldh+k-1]) + } + if k >= ktop+2 { + tst1 += math.Abs(h[k*ldh+k-2]) + } + if k >= ktop+3 { + tst1 += math.Abs(h[k*ldh+k-3]) + } + if k <= kbot-2 { + tst1 += math.Abs(h[(k+2)*ldh+k+1]) + } + if k <= kbot-3 { + tst1 += math.Abs(h[(k+3)*ldh+k+1]) + } + if k <= kbot-4 { + tst1 += math.Abs(h[(k+4)*ldh+k+1]) + } } - if k <= kbot-2 { - tst1 += math.Abs(h[(k+2)*ldh+k+1]) + if math.Abs(h[(k+1)*ldh+k]) <= math.Max(smlnum, ulp*tst1) { + h12 := math.Max(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) + h21 := math.Min(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) + h11 := math.Max(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) + h22 := math.Min(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) + scl := h11 + h12 + tst2 := h22 * (h11 / scl) + if tst2 == 0 || h21*(h12/scl) <= math.Max(smlnum, ulp*tst2) { + h[(k+1)*ldh+k] = 0 + } } - if k <= kbot-3 { - tst1 += math.Abs(h[(k+3)*ldh+k+1]) - } - if k <= kbot-4 { - tst1 += math.Abs(h[(k+4)*ldh+k+1]) + } + } + // Multiply H by reflections from the left. + var jbot int + switch { + case accum: + jbot = min(ndcol, kbot) + case wantt: + jbot = n - 1 + default: + jbot = kbot + } + for m := mbot; m >= mtop; m-- { + k := krcol + 2*m + t1 := v[m*ldv] + t2 := t1 * v[m*ldv+1] + t3 := t1 * v[m*ldv+2] + for j := max(ktop, krcol+2*(m+1)); j <= jbot; j++ { + refsum := h[(k+1)*ldh+j] + v[m*ldv+1]*h[(k+2)*ldh+j] + v[m*ldv+2]*h[(k+3)*ldh+j] + h[(k+1)*ldh+j] -= refsum * t1 + h[(k+2)*ldh+j] -= refsum * t2 + h[(k+3)*ldh+j] -= refsum * t3 + } + } + // Accumulate orthogonal transformations. + if accum { + // Accumulate U. If necessary, update Z later with an + // efficient matrix-matrix multiply. + for m := mbot; m >= mtop; m-- { + k := krcol + 2*m + kms := k - incol - 1 + i2 := max(0, ktop-incol-1) + i2 = max(i2, kms-(krcol-incol)) + i4 := min(kdu, krcol+2*mbot-incol+5) + t1 := v[m*ldv] + t2 := t1 * v[m*ldv+1] + t3 := t1 * v[m*ldv+2] + for j := i2; j < i4; j++ { + refsum := u[j*ldu+kms+1] + v[m*ldv+1]*u[j*ldu+kms+2] + v[m*ldv+2]*u[j*ldu+kms+3] + u[j*ldu+kms+1] -= refsum * t1 + u[j*ldu+kms+2] -= refsum * t2 + u[j*ldu+kms+3] -= refsum * t3 } } - if math.Abs(h[(k+1)*ldh+k]) <= math.Max(smlnum, ulp*tst1) { - h12 := math.Max(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) - h21 := math.Min(math.Abs(h[(k+1)*ldh+k]), math.Abs(h[k*ldh+k+1])) - h11 := math.Max(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) - h22 := math.Min(math.Abs(h[(k+1)*ldh+k+1]), math.Abs(h[k*ldh+k]-h[(k+1)*ldh+k+1])) - scl := h11 + h12 - tst2 := h22 * (h11 / scl) - if tst2 == 0 || h21*(h12/scl) <= math.Max(smlnum, ulp*tst2) { - h[(k+1)*ldh+k] = 0 + } else if wantz { + // U is not accumulated, so update Z now by multiplying by + // reflections from the right. + for m := mbot; m >= mtop; m-- { + k := krcol + 2*m + t1 := v[m*ldv] + t2 := t1 * v[m*ldv+1] + t3 := t1 * v[m*ldv+2] + for j := iloz; j <= ihiz; j++ { + refsum := z[j*ldz+k+1] + v[m*ldv+1]*z[j*ldz+k+2] + v[m*ldv+2]*z[j*ldz+k+3] + z[j*ldz+k+1] -= refsum * t1 + z[j*ldz+k+2] -= refsum * t2 + z[j*ldz+k+3] -= refsum * t3 } } } - - // Fill in the last row of each bulge. - mend = min(nbmps, (kbot-krcol-1)/3) - 1 - for m := mtop; m <= mend; m++ { - k := krcol + 3*m - refsum := v[m*ldv] * v[m*ldv+2] * h[(k+4)*ldh+k+3] - h[(k+4)*ldh+k+1] = -refsum - h[(k+4)*ldh+k+2] = -refsum * v[m*ldv+1] - h[(k+4)*ldh+k+3] -= refsum * v[m*ldv+2] - } } - // Use U (if accumulated) to update far-from-diagonal entries in H. // If required, use U to update Z as well. if !accum { continue } - var jtop, jbot int + jtop, jbot := ktop, kbot if wantt { jtop = 0 jbot = n - 1 - } else { - jtop = ktop - jbot = kbot } bi := blas64.Implementation() - if !blk22 || incol < ktop || kbot < ndcol || ns <= 2 { - // Updates not exploiting the 2×2 block structure of U. k0 and nu keep track - // of the location and size of U in the special cases of introducing bulges - // and chasing bulges off the bottom. In these special cases and in case the - // number of shifts is ns = 2, there is no 2×2 block structure to exploit. - - k0 := max(0, ktop-incol-1) - nu := kdu - max(0, ndcol-kbot) - k0 - - // Horizontal multiply. - for jcol := min(ndcol, kbot) + 1; jcol <= jbot; jcol += nh { - jlen := min(nh, jbot-jcol+1) - bi.Dgemm(blas.Trans, blas.NoTrans, nu, jlen, nu, - 1, u[k0*ldu+k0:], ldu, - h[(incol+k0+1)*ldh+jcol:], ldh, - 0, wh, ldwh) - impl.Dlacpy(blas.All, nu, jlen, wh, ldwh, h[(incol+k0+1)*ldh+jcol:], ldh) - } - - // Vertical multiply. - for jrow := jtop; jrow <= max(ktop, incol)-1; jrow += nv { - jlen := min(nv, max(ktop, incol)-jrow) - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, nu, nu, - 1, h[jrow*ldh+incol+k0+1:], ldh, - u[k0*ldu+k0:], ldu, - 0, wv, ldwv) - impl.Dlacpy(blas.All, jlen, nu, wv, ldwv, h[jrow*ldh+incol+k0+1:], ldh) - } - - // Z multiply (also vertical). - if wantz { - for jrow := iloz; jrow <= ihiz; jrow += nv { - jlen := min(nv, ihiz-jrow+1) - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, nu, nu, - 1, z[jrow*ldz+incol+k0+1:], ldz, - u[k0*ldu+k0:], ldu, - 0, wv, ldwv) - impl.Dlacpy(blas.All, jlen, nu, wv, ldwv, z[jrow*ldz+incol+k0+1:], ldz) - } - } - - continue - } - - // Updates exploiting U's 2×2 block structure. - - // i2, i4, j2, j4 are the last rows and columns of the blocks. - i2 := (kdu + 1) / 2 - i4 := kdu - j2 := i4 - i2 - j4 := kdu - - // kzs and knz deal with the band of zeros along the diagonal of one of the - // triangular blocks. - kzs := (j4 - j2) - (ns + 1) - knz := ns + 1 - + k1 := max(0, ktop-incol-1) + nu := kdu - max(0, ndcol-kbot) - k1 // Horizontal multiply. for jcol := min(ndcol, kbot) + 1; jcol <= jbot; jcol += nh { jlen := min(nh, jbot-jcol+1) - - // Copy bottom of H to top+kzs of scratch (the first kzs - // rows get multiplied by zero). - impl.Dlacpy(blas.All, knz, jlen, h[(incol+1+j2)*ldh+jcol:], ldh, wh[kzs*ldwh:], ldwh) - - // Multiply by U21ᵀ. - impl.Dlaset(blas.All, kzs, jlen, 0, 0, wh, ldwh) - bi.Dtrmm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit, knz, jlen, - 1, u[j2*ldu+kzs:], ldu, wh[kzs*ldwh:], ldwh) - - // Multiply top of H by U11ᵀ. - bi.Dgemm(blas.Trans, blas.NoTrans, i2, jlen, j2, - 1, u, ldu, h[(incol+1)*ldh+jcol:], ldh, - 1, wh, ldwh) - - // Copy top of H to bottom of WH. - impl.Dlacpy(blas.All, j2, jlen, h[(incol+1)*ldh+jcol:], ldh, wh[i2*ldwh:], ldwh) - - // Multiply by U21ᵀ. - bi.Dtrmm(blas.Left, blas.Lower, blas.Trans, blas.NonUnit, j2, jlen, - 1, u[i2:], ldu, wh[i2*ldwh:], ldwh) - - // Multiply by U22. - bi.Dgemm(blas.Trans, blas.NoTrans, i4-i2, jlen, j4-j2, - 1, u[j2*ldu+i2:], ldu, h[(incol+1+j2)*ldh+jcol:], ldh, - 1, wh[i2*ldwh:], ldwh) - - // Copy it back. - impl.Dlacpy(blas.All, kdu, jlen, wh, ldwh, h[(incol+1)*ldh+jcol:], ldh) + bi.Dgemm(blas.Trans, blas.NoTrans, nu, jlen, nu, + 1, u[k1*ldu+k1:], ldu, + h[(incol+k1+1)*ldh+jcol:], ldh, + 0, wh, ldwh) + impl.Dlacpy(blas.All, nu, jlen, wh, ldwh, h[(incol+k1+1)*ldh+jcol:], ldh) } - // Vertical multiply. - for jrow := jtop; jrow <= max(incol, ktop)-1; jrow += nv { - jlen := min(nv, max(incol, ktop)-jrow) - - // Copy right of H to scratch (the first kzs columns get multiplied - // by zero). - impl.Dlacpy(blas.All, jlen, knz, h[jrow*ldh+incol+1+j2:], ldh, wv[kzs:], ldwv) - - // Multiply by U21. - impl.Dlaset(blas.All, jlen, kzs, 0, 0, wv, ldwv) - bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.NonUnit, jlen, knz, - 1, u[j2*ldu+kzs:], ldu, wv[kzs:], ldwv) - - // Multiply by U11. - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, i2, j2, - 1, h[jrow*ldh+incol+1:], ldh, u, ldu, - 1, wv, ldwv) - - // Copy left of H to right of scratch. - impl.Dlacpy(blas.All, jlen, j2, h[jrow*ldh+incol+1:], ldh, wv[i2:], ldwv) - - // Multiply by U21. - bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.NonUnit, jlen, i4-i2, - 1, u[i2:], ldu, wv[i2:], ldwv) - - // Multiply by U22. - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, i4-i2, j4-j2, - 1, h[jrow*ldh+incol+1+j2:], ldh, u[j2*ldu+i2:], ldu, - 1, wv[i2:], ldwv) - - // Copy it back. - impl.Dlacpy(blas.All, jlen, kdu, wv, ldwv, h[jrow*ldh+incol+1:], ldh) + for jrow := jtop; jrow < max(ktop, incol); jrow += nv { + jlen := min(nv, max(ktop, incol)-jrow) + bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, nu, nu, + 1, h[jrow*ldh+incol+k1+1:], ldh, + u[k1*ldu+k1:], ldu, + 0, wv, ldwv) + impl.Dlacpy(blas.All, jlen, nu, wv, ldwv, h[jrow*ldh+incol+k1+1:], ldh) } - - if !wantz { - continue - } - // Multiply Z (also vertical). - for jrow := iloz; jrow <= ihiz; jrow += nv { - jlen := min(nv, ihiz-jrow+1) - - // Copy right of Z to left of scratch (first kzs columns get - // multiplied by zero). - impl.Dlacpy(blas.All, jlen, knz, z[jrow*ldz+incol+1+j2:], ldz, wv[kzs:], ldwv) - - // Multiply by U12. - impl.Dlaset(blas.All, jlen, kzs, 0, 0, wv, ldwv) - bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.NonUnit, jlen, knz, - 1, u[j2*ldu+kzs:], ldu, wv[kzs:], ldwv) - - // Multiply by U11. - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, i2, j2, - 1, z[jrow*ldz+incol+1:], ldz, u, ldu, - 1, wv, ldwv) - - // Copy left of Z to right of scratch. - impl.Dlacpy(blas.All, jlen, j2, z[jrow*ldz+incol+1:], ldz, wv[i2:], ldwv) - - // Multiply by U21. - bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.NonUnit, jlen, i4-i2, - 1, u[i2:], ldu, wv[i2:], ldwv) - - // Multiply by U22. - bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, i4-i2, j4-j2, - 1, z[jrow*ldz+incol+1+j2:], ldz, u[j2*ldu+i2:], ldu, - 1, wv[i2:], ldwv) - - // Copy the result back to Z. - impl.Dlacpy(blas.All, jlen, kdu, wv, ldwv, z[jrow*ldz+incol+1:], ldz) + // Z multiply (also vertical). + if wantz { + for jrow := iloz; jrow <= ihiz; jrow += nv { + jlen := min(nv, ihiz-jrow+1) + bi.Dgemm(blas.NoTrans, blas.NoTrans, jlen, nu, nu, + 1, z[jrow*ldz+incol+k1+1:], ldz, + u[k1*ldu+k1:], ldu, + 0, wv, ldwv) + impl.Dlacpy(blas.All, jlen, nu, wv, ldwv, z[jrow*ldz+incol+k1+1:], ldz) + } } } } diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlassq.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlassq.go index 3d982c3c..8f8cda87 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlassq.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlassq.go @@ -77,21 +77,25 @@ func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, s ax := scale * math.Sqrt(sumsq) switch { case ax > dtbig: - // We assume scale >= sqrt( TINY*EPS ) / dsbig, that is, if the - // scaled sum is big then its scaling factor should not be too - // small. - v := scale * dsbig - abig += (v * v) * sumsq + if scale > 1 { + scale *= dsbig + abig += scale * (scale * sumsq) + } else { + // sumsq > dtbig^2 => (dsbig * (dsbig * sumsq)) is representable. + abig += scale * (scale * (dsbig * (dsbig * sumsq))) + } case ax < dtsml: if !isBig { - // We assume scale <= sqrt( HUGE ) / dssml, that is, if the - // scaled sum is small then its scaling factor should not be too - // big. - v := scale * dssml - asml += (v * v) * sumsq + if scale < 1 { + scale *= dssml + asml += scale * (scale * sumsq) + } else { + // sumsq < dtsml^2 => (dssml * (dssml * sumsq)) is representable. + asml += scale * (scale * (dssml * (dssml * sumsq))) + } } default: - amed += ax * ax + amed += scale * (scale * sumsq) } } // Combine abig and amed or amed and asml if more than one accumulator was diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatbs.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlatbs.go index 19300faf..e0e809cf 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatbs.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlatbs.go @@ -52,7 +52,7 @@ func (Implementation) Dlatbs(uplo blas.Uplo, trans blas.Transpose, diag blas.Dia // Quick return if possible. if n == 0 { - return 0 + return 1 } switch { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go index 37ac2fe7..f13b7d57 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dlatrs.go @@ -9,6 +9,7 @@ import ( "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" + "gonum.org/v1/gonum/lapack" ) // Dlatrs solves a triangular system of equations scaled to prevent overflow. It @@ -45,7 +46,7 @@ func (impl Implementation) Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag bla // Quick return if possible. if n == 0 { - return 0 + return 1 } switch { @@ -81,13 +82,61 @@ func (impl Implementation) Dlatrs(uplo blas.Uplo, trans blas.Transpose, diag bla } // Scale the column norms by tscal if the maximum element in cnorm is greater than bignum. imax := bi.Idamax(n, cnorm, 1) - tmax := cnorm[imax] var tscal float64 - if tmax <= bignum { + if cnorm[imax] <= bignum { tscal = 1 } else { - tscal = 1 / (smlnum * tmax) - bi.Dscal(n, tscal, cnorm, 1) + tmax := cnorm[imax] + // Avoid NaN generation if entries in cnorm exceed the overflow + // threshold. + if tmax <= math.MaxFloat64 { + // Case 1: All entries in cnorm are valid floating-point numbers. + tscal = 1 / (smlnum * tmax) + bi.Dscal(n, tscal, cnorm, 1) + } else { + // Case 2: At least one column norm of A cannot be represented as + // floating-point number. Find the offdiagonal entry A[i,j] with the + // largest absolute value. If this entry is not +/- Infinity, use + // this value as tscal. + tmax = 0 + if upper { + // A is upper triangular. + for j := 1; j < n; j++ { + tmax = math.Max(impl.Dlange(lapack.MaxAbs, j, 1, a[j:], lda, nil), tmax) + } + } else { + // A is lower triangular. + for j := 0; j < n-1; j++ { + tmax = math.Max(impl.Dlange(lapack.MaxAbs, n-j-1, 1, a[(j+1)*lda+j:], lda, nil), tmax) + } + } + if tmax <= math.MaxFloat64 { + tscal = 1 / (smlnum * tmax) + for j := 0; j < n; j++ { + if cnorm[j] <= math.MaxFloat64 { + cnorm[j] *= tscal + } else { + // Recompute the 1-norm without introducing Infinity in + // the summation. + cnorm[j] = 0 + if upper { + for i := 0; i < j; i++ { + cnorm[j] += tscal * math.Abs(a[i*lda+j]) + } + } else { + for i := j + 1; i < n; i++ { + cnorm[j] += tscal * math.Abs(a[i*lda+j]) + } + } + } + } + } else { + // At least one entry of A is not a valid floating-point entry. + // Rely on Dtrsv to propagate Inf and NaN. + bi.Dtrsv(uplo, trans, diag, n, a, lda, x, 1) + return + } + } } // Compute a bound on the computed solution vector to see if bi.Dtrsv can be used. diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorg2r.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorg2r.go index 36addec7..c56f24cb 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorg2r.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorg2r.go @@ -14,7 +14,7 @@ import ( // // Q = H_0 * H_1 * ... * H_{k-1} // -// len(tau) >= k, 0 <= k <= n, 0 <= n <= m, len(work) >= n. +// len(tau) = k, 0 <= k <= n, 0 <= n <= m, len(work) >= n. // Dorg2r will panic if these conditions are not met. // // Dorg2r is an internal routine. It is exported for testing purposes. @@ -41,8 +41,8 @@ func (impl Implementation) Dorg2r(m, n, k int, a []float64, lda int, tau []float switch { case len(a) < (m-1)*lda+n: panic(shortA) - case len(tau) < k: - panic(shortTau) + case len(tau) != k: + panic(badLenTau) case len(work) < n: panic(shortWork) } diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgbr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgbr.go index 8a4fe2b5..35535100 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgbr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgbr.go @@ -91,7 +91,7 @@ func (impl Implementation) Dorgbr(vect lapack.GenOrtho, m, n, k int, a []float64 if wantq { // Form Q, determined by a call to Dgebrd to reduce an m×k matrix. if m >= k { - impl.Dorgqr(m, n, k, a, lda, tau, work, lwork) + impl.Dorgqr(m, n, k, a, lda, tau[:k], work, lwork) } else { // Shift the vectors which define the elementary reflectors one // column to the right, and set the first row and column of Q to @@ -108,7 +108,7 @@ func (impl Implementation) Dorgbr(vect lapack.GenOrtho, m, n, k int, a []float64 } if m > 1 { // Form Q[1:m-1, 1:m-1] - impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork) + impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau[:m-1], work, lwork) } } } else { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgl2.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgl2.go index ea6dbe52..6dd9a888 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgl2.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgl2.go @@ -9,13 +9,19 @@ import ( "gonum.org/v1/gonum/blas/blas64" ) -// Dorgl2 generates an m×n matrix Q with orthonormal rows defined by the -// first m rows product of elementary reflectors as computed by Dgelqf. +// Dorgl2 generates an m×n matrix Q with orthonormal rows defined as the first m +// rows of a product of k elementary reflectors of order n // -// Q = H_0 * H_1 * ... * H_{k-1} +// Q = H_{k-1} * ... * H_0 // -// len(tau) >= k, 0 <= k <= m, 0 <= m <= n, len(work) >= m. -// Dorgl2 will panic if these conditions are not met. +// as returned by Dgelqf. +// +// On entry, tau and the first k rows of A must contain the scalar factors and +// the vectors, respectively, which define the elementary reflectors H_i, +// i=0,...,k-1, as returned by Dgelqf. On return, A contains the matrix Q. +// +// tau must have length at least k, work must have length at least m, and it +// must hold that 0 <= k <= m <= n, otherwise Dorgl2 will panic. // // Dorgl2 is an internal routine. It is exported for testing purposes. func (impl Implementation) Dorgl2(m, n, k int, a []float64, lda int, tau, work []float64) { @@ -28,7 +34,7 @@ func (impl Implementation) Dorgl2(m, n, k int, a []float64, lda int, tau, work [ panic(kLT0) case k > m: panic(kGTM) - case lda < max(1, m): + case lda < max(1, n): panic(badLdA) } diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorglq.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorglq.go index 1128eb30..d6b3aadf 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorglq.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorglq.go @@ -9,24 +9,24 @@ import ( "gonum.org/v1/gonum/lapack" ) -// Dorglq generates an m×n matrix Q with orthonormal columns defined by the -// product of elementary reflectors as computed by Dgelqf. +// Dorglq generates an m×n matrix Q with orthonormal rows defined as the first m +// rows of a product of k elementary reflectors of order n // -// Q = H_0 * H_1 * ... * H_{k-1} +// Q = H_{k-1} * ... * H_0 // -// Dorglq is the blocked version of Dorgl2 that makes greater use of level-3 BLAS -// routines. +// as returned by Dgelqf. // -// len(tau) >= k, 0 <= k <= m, and 0 <= m <= n. +// On entry, tau and the first k rows of A must contain the scalar factors and +// the vectors, respectively, which define the elementary reflectors H_i, +// i=0,...,k-1, as returned by Dgelqf. On return, A contains the matrix Q. // -// work is temporary storage, and lwork specifies the usable memory length. At minimum, -// lwork >= m, and the amount of blocking is limited by the usable length. -// If lwork == -1, instead of computing Dorglq the optimal work length is stored -// into work[0]. +// tau must have length at least k, work must have length at least lwork and +// lwork must be at least max(1,m). On return, optimal value of lwork will be +// stored in work[0]. It must also hold that 0 <= k <= m <= n, otherwise Dorglq +// will panic. // -// Dorglq will panic if the conditions on input values are not met. -// -// Dorglq is an internal routine. It is exported for testing purposes. +// If lwork == -1, instead of performing Dorglq, the function only calculates +// the optimal value of lwork and stores it into work[0]. func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) { switch { case m < 0: diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgqr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgqr.go index db1a4807..a1e0fa87 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgqr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgqr.go @@ -18,7 +18,7 @@ import ( // Dorgqr is the blocked version of Dorg2r that makes greater use of level-3 BLAS // routines. // -// The length of tau must be at least k, and the length of work must be at least n. +// The length of tau must be k, and the length of work must be at least n. // It also must be that 0 <= k <= n and 0 <= n <= m. // // work is temporary storage, and lwork specifies the usable memory length. At @@ -70,8 +70,8 @@ func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work [ switch { case len(a) < (m-1)*lda+n: panic(shortA) - case len(tau) < k: - panic(shortTau) + case len(tau) != k: + panic(badLenTau) } nbmin := 2 // Minimum block size @@ -123,7 +123,7 @@ func (impl Implementation) Dorgqr(m, n, k int, a []float64, lda int, tau, work [ a[i*lda+i+ib:], lda, work[ib*ldwork:], ldwork) } - impl.Dorg2r(m-i, ib, ib, a[i*lda+i:], lda, tau[i:], work) + impl.Dorg2r(m-i, ib, ib, a[i*lda+i:], lda, tau[i:i+ib], work) // Set rows 0:i-1 of current block to zero. for j := i; j < i+ib; j++ { for l := 0; l < i; l++ { diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgtr.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgtr.go index d2cef99a..7021ae53 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorgtr.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorgtr.go @@ -99,7 +99,7 @@ func (impl Implementation) Dorgtr(uplo blas.Uplo, n int, a []float64, lda int, t } if n > 1 { // Generate Q[1:n, 1:n]. - impl.Dorgqr(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork) + impl.Dorgqr(n-1, n-1, n-1, a[lda+1:], lda, tau[:n-1], work, lwork) } } work[0] = float64(lworkopt) diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dorm2r.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dorm2r.go index e9293131..aea77a70 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/dorm2r.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dorm2r.go @@ -17,7 +17,7 @@ import "gonum.org/v1/gonum/blas" // If side == blas.Left, a is a matrix of size m×k, and if side == blas.Right // a is of size n×k. // -// tau contains the Householder factors and is of length at least k and this function +// tau contains the Householder factors and must have length k and this function // will panic otherwise. // // work is temporary storage of length at least n if side == blas.Left @@ -59,8 +59,8 @@ func (impl Implementation) Dorm2r(side blas.Side, trans blas.Transpose, m, n, k panic(shortA) case len(c) < (m-1)*ldc+n: panic(shortC) - case len(tau) < k: - panic(shortTau) + case len(tau) != k: + panic(badLenTau) case left && len(work) < n: panic(shortWork) case !left && len(work) < m: diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dptcon.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dptcon.go new file mode 100644 index 00000000..cd41e317 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dptcon.go @@ -0,0 +1,99 @@ +// Copyright ©2023 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/blas64" +) + +// Dptcon computes and returns the reciprocal of the condition number (in the +// 1-norm) of a symmetric positive definite tridiagonal matrix A using the +// factorization A = L*D*Lᵀ or A = Uᵀ*D*U computed by Dpttrf. +// +// The reciprocal of the condition number is computed as +// +// rcond = 1 / (anorm * ‖A⁻¹‖) +// +// and ‖A⁻¹‖ is computed by a direct method. +// +// d and e contain, respectively, the n diagonal elements of the diagonal matrix +// D and the (n-1) off-diagonal elements of the unit bidiagonal factor U or L +// from the factorization of A, as computed by Dpttrf. +// +// anorm is the 1-norm of the original matrix A. +// +// work must have length n, otherwise Dptcon will panic. +func (impl Implementation) Dptcon(n int, d, e []float64, anorm float64, work []float64) (rcond float64) { + switch { + case n < 0: + panic(nLT0) + case anorm < 0: + panic(badNorm) + } + + // Quick return if possible. + if n == 0 { + return 1 + } + + switch { + case len(d) < n: + panic(shortD) + case len(e) < n-1: + panic(shortE) + case len(work) < n: + panic(shortWork) + } + + // Quick return if possible. + switch { + case anorm == 0: + return 0 + case math.IsNaN(anorm): + // Propagate NaN. + return anorm + case math.IsInf(anorm, 1): + return 0 + } + + // Check that d[0:n] is positive. + for _, di := range d[:n] { + if di <= 0 { + return 0 + } + } + + // Solve M(A) * x = e, where M(A) = (m[i,j]) is given by + // + // m[i,j] = abs(A[i,j]), i == j, + // m[i,j] = -abs(A[i,j]), i != j, + // + // and e = [1,1,...,1]ᵀ. Note M(A) = M(L)*D*M(L)ᵀ. + + // Solve M(L) * b = e. + work[0] = 1 + for i := 1; i < n; i++ { + work[i] = 1 + work[i-1]*math.Abs(e[i-1]) + } + + // Solve D * M(L)ᵀ * x = b. + work[n-1] /= d[n-1] + for i := n - 2; i >= 0; i-- { + work[i] = work[i]/d[i] + work[i+1]*math.Abs(e[i]) + } + + // Compute ainvnm = max(x[i]), 0<=i 0 +} diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dpttrs.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dpttrs.go new file mode 100644 index 00000000..7bdee6f9 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dpttrs.go @@ -0,0 +1,51 @@ +// Copyright ©2023 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 + +// Dpttrs solves a tridiagonal system of the form +// +// A * X = B +// +// using the L*D*Lᵀ factorization of A computed by Dpttrf. D is a diagonal +// matrix specified in d, L is a unit bidiagonal matrix whose subdiagonal is +// specified in e, and X and B are n×nrhs matrices. +func (impl Implementation) Dpttrs(n, nrhs int, d, e []float64, b []float64, ldb int) { + switch { + case n < 0: + panic(nLT0) + case nrhs < 0: + panic(nrhsLT0) + case ldb < max(1, nrhs): + panic(badLdB) + } + + // Quick return if possible. + if n == 0 || nrhs == 0 { + return + } + + switch { + case len(d) < n: + panic(shortD) + case len(e) < n-1: + panic(shortE) + case len(b) < (n-1)*ldb+nrhs: + panic(shortB) + } + + nb := 1 + if nrhs > 1 { + nb = max(1, impl.Ilaenv(1, "DPTTRS", " ", n, nrhs, -1, -1)) + } + + if nb >= nrhs { + impl.dptts2(n, nrhs, d, e, b, ldb) + } else { + for j := 0; j < nrhs; j += nb { + jb := min(nrhs-j, nb) + impl.dptts2(n, jb, d, e, b[j:], ldb) + } + } +} diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/dptts2.go b/vendor/gonum.org/v1/gonum/lapack/gonum/dptts2.go new file mode 100644 index 00000000..ff1df168 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/dptts2.go @@ -0,0 +1,39 @@ +// Copyright ©2023 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 "gonum.org/v1/gonum/blas/blas64" + +// dptts2 solves a tridiagonal system of the form +// +// A * X = B +// +// using the L*D*Lᵀ factorization of A computed by Dpttrf. D is a diagonal +// matrix specified in d, L is a unit bidiagonal matrix whose subdiagonal is +// specified in e, and X and B are n×nrhs matrices. +func (impl Implementation) dptts2(n, nrhs int, d, e []float64, b []float64, ldb int) { + // Quick return if possible. + if n <= 1 { + if n == 1 { + bi := blas64.Implementation() + bi.Dscal(nrhs, 1/d[0], b, 1) + } + return + } + + // Solve A * X = B using the factorization A = L*D*Lᵀ, overwriting each + // right hand side vector with its solution. + for j := 0; j < nrhs; j++ { + // Solve L * x = b. + for i := 1; i < n; i++ { + b[i*ldb+j] -= b[(i-1)*ldb+j] * e[i-1] + } + // Solve D * Lᵀ * x = b. + b[(n-1)*ldb+j] /= d[n-1] + for i := n - 2; i >= 0; i-- { + b[i*ldb+j] = b[i*ldb+j]/d[i] - b[(i+1)*ldb+j]*e[i] + } + } +} diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/errors.go b/vendor/gonum.org/v1/gonum/lapack/gonum/errors.go index 47652841..711cc2d5 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/errors.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/errors.go @@ -21,6 +21,7 @@ const ( badMatrixType = "lapack: bad MatrixType" badMaximizeNormXJob = "lapack: bad MaximizeNormXJob" badNorm = "lapack: bad Norm" + badOrthoComp = "lapack: bad OrthoComp" badPivot = "lapack: bad Pivot" badRightEVJob = "lapack: bad RightEVJob" badSVDJob = "lapack: bad SVDJob" diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/ilaenv.go b/vendor/gonum.org/v1/gonum/lapack/gonum/ilaenv.go index c79597e3..fc70806c 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/ilaenv.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/ilaenv.go @@ -176,6 +176,13 @@ func (impl Implementation) Ilaenv(ispec int, name string, opts string, n1, n2, n } return 32 } + case "PT": + switch c3 { + default: + panic(badName) + case "TRS": + return 1 + } case "TR": switch c3 { default: diff --git a/vendor/gonum.org/v1/gonum/lapack/gonum/lapack.go b/vendor/gonum.org/v1/gonum/lapack/gonum/lapack.go index fef4f558..5daefc58 100644 --- a/vendor/gonum.org/v1/gonum/lapack/gonum/lapack.go +++ b/vendor/gonum.org/v1/gonum/lapack/gonum/lapack.go @@ -13,20 +13,6 @@ type Implementation struct{} var _ lapack.Float64 = Implementation{} -func min(a, b int) int { - if a < b { - return a - } - return b -} - -func max(a, b int) int { - if a > b { - return a - } - return b -} - func abs(a int) int { if a < 0 { return -a diff --git a/vendor/gonum.org/v1/gonum/lapack/lapack.go b/vendor/gonum.org/v1/gonum/lapack/lapack.go index 5f60438f..60ef1c24 100644 --- a/vendor/gonum.org/v1/gonum/lapack/lapack.go +++ b/vendor/gonum.org/v1/gonum/lapack/lapack.go @@ -15,6 +15,7 @@ type Float64 interface { Dgeev(jobvl LeftEVJob, jobvr RightEVJob, n int, a []float64, lda int, wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) (first int) Dgels(trans blas.Transpose, m, n, nrhs int, a []float64, lda int, b []float64, ldb int, work []float64, lwork int) bool Dgelqf(m, n int, a []float64, lda int, tau, work []float64, lwork int) + Dgeqp3(m, n int, a []float64, lda int, jpvt []int, tau, work []float64, lwork int) Dgeqrf(m, n int, a []float64, lda int, tau, work []float64, lwork int) Dgesvd(jobU, jobVT SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool) Dgetrf(m, n int, a []float64, lda int, ipiv []int) (ok bool) @@ -26,7 +27,9 @@ type Float64 interface { Dlansy(norm MatrixNorm, uplo blas.Uplo, n int, a []float64, lda int, work []float64) float64 Dlapmr(forward bool, m, n int, x []float64, ldx int, k []int) Dlapmt(forward bool, m, n int, x []float64, ldx int, k []int) + Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) + Dorglq(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) float64 Dpbtrf(uplo blas.Uplo, n, kd int, ab []float64, ldab int) (ok bool) @@ -218,7 +221,7 @@ const ( EVSelected EVHowMany = 'S' // Compute selected right and/or left eigenvectors. ) -// MaximizeNormX specifies the heuristic method for computing a contribution to +// MaximizeNormXJob specifies the heuristic method for computing a contribution to // the reciprocal Dif-estimate in Dlatdf. type MaximizeNormXJob byte @@ -226,3 +229,12 @@ const ( LocalLookAhead MaximizeNormXJob = 0 // Solve Z*x=h-f where h is a vector of ±1. NormalizedNullVector MaximizeNormXJob = 2 // Compute an approximate null-vector e of Z, normalize e and solve Z*x=±e-f. ) + +// OrthoComp specifies whether and how the orthogonal matrix is computed in Dgghrd. +type OrthoComp byte + +const ( + OrthoNone OrthoComp = 'N' // Do not compute the orthogonal matrix. + OrthoExplicit OrthoComp = 'I' // The orthogonal matrix is formed explicitly and returned in the argument. + OrthoPostmul OrthoComp = 'V' // The orthogonal matrix is post-multiplied into the matrix stored in the argument on entry. +) diff --git a/vendor/gonum.org/v1/gonum/lapack/lapack64/lapack64.go b/vendor/gonum.org/v1/gonum/lapack/lapack64/lapack64.go index acb62da4..1b4c1734 100644 --- a/vendor/gonum.org/v1/gonum/lapack/lapack64/lapack64.go +++ b/vendor/gonum.org/v1/gonum/lapack/lapack64/lapack64.go @@ -27,13 +27,6 @@ type Tridiagonal struct { DU []float64 } -func max(a, b int) int { - if a > b { - return a - } - return b -} - // Potrf computes the Cholesky factorization of a. // The factorization has the form // @@ -222,11 +215,52 @@ func Gels(trans blas.Transpose, a blas64.General, b blas64.General, work []float return lapack64.Dgels(trans, a.Rows, a.Cols, b.Cols, a.Data, max(1, a.Stride), b.Data, max(1, b.Stride), work, lwork) } +// Geqp3 computes a QR factorization with column pivoting of the m×n matrix A: +// +// A*P = Q*R +// +// where P is a permutation matrix, Q is an orthogonal matrix and R is a +// min(m,n)×n upper trapezoidal matrix. +// +// On return, the upper triangle of A contains the matrix R. The elements below +// the diagonal together with tau represent the matrix Q as a product of +// elementary reflectors +// +// Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n). +// +// Each H_i has the form +// +// H_i = I - tau * v * vᵀ +// +// where tau is a scalar and v is a vector with v[0:i] = 0 and v[i] = 1; +// v[i+1:m] is stored on exit in A[i+1:m,i], and tau in tau[i]. +// +// jpvt specifies a column pivot to be applied to A. On entry, if jpvt[j] is at +// least zero, the jth column of A is permuted to the front of A*P (a leading +// column), if jpvt[j] is -1 the jth column of A is a free column. If jpvt[j] < +// -1, Geqp3 will panic. On return, jpvt holds the permutation that was applied; +// the jth column of A*P was the jpvt[j] column of A. jpvt must have length n or +// Geqp3 will panic. +// +// tau holds the scalar factors of the elementary reflectors. It must have +// length min(m,n), otherwise Geqp3 will panic. +// +// work must have length at least max(1,lwork), and lwork must be at least +// 3*n+1, otherwise Geqp3 will panic. For optimal performance lwork must be at +// least 2*n+(n+1)*nb, where nb is the optimal blocksize. On return, work[0] +// will contain the optimal value of lwork. +// +// If lwork == -1, instead of performing Geqp3, only the optimal value of lwork +// will be stored in work[0]. +func Geqp3(a blas64.General, jpvt []int, tau, work []float64, lwork int) { + lapack64.Dgeqp3(a.Rows, a.Cols, a.Data, max(1, a.Stride), jpvt, tau, work, lwork) +} + // Geqrf computes the QR factorization of the m×n matrix A using a blocked // algorithm. A is modified to contain the information to construct Q and R. // The upper triangle of a contains the matrix R. The lower triangular elements // (not including the diagonal) contain the elementary reflectors. tau is modified -// to contain the reflector scales. tau must have length at least min(m,n), and +// to contain the reflector scales. tau must have length min(m,n), and // this function will panic otherwise. // // The ith elementary reflector can be explicitly constructed by first extracting @@ -238,7 +272,7 @@ func Gels(trans blas.Transpose, a blas64.General, b blas64.General, work []float // // and computing H_i = I - tau[i] * v * vᵀ. // -// The orthonormal matrix Q can be constucted from a product of these elementary +// The orthonormal matrix Q can be constructed from a product of these elementary // reflectors, Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n). // // Work is temporary storage, and lwork specifies the usable memory length. @@ -320,25 +354,27 @@ func Gesvd(jobU, jobVT lapack.SVDJob, a, u, vt blas64.General, s, work []float64 return lapack64.Dgesvd(jobU, jobVT, a.Rows, a.Cols, a.Data, max(1, a.Stride), s, u.Data, max(1, u.Stride), vt.Data, max(1, vt.Stride), work, lwork) } -// Getrf computes the LU decomposition of the m×n matrix A. +// Getrf computes the LU decomposition of an m×n matrix A using partial +// pivoting with row interchanges. +// // The LU decomposition is a factorization of A into // // A = P * L * U // -// where P is a permutation matrix, L is a unit lower triangular matrix, and -// U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored -// in place into a. +// where P is a permutation matrix, L is a lower triangular with unit diagonal +// elements (lower trapezoidal if m > n), and U is upper triangular (upper +// trapezoidal if m < n). // -// ipiv is a permutation vector. It indicates that row i of the matrix was -// changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic -// otherwise. ipiv is zero-indexed. +// On entry, a contains the matrix A. On return, L and U are stored in place +// into a, and P is represented by ipiv. // -// Getrf is the blocked version of the algorithm. +// ipiv contains a sequence of row swaps. It indicates that row i of the matrix +// was interchanged with ipiv[i]. ipiv must have length min(m,n), and Getrf will +// panic otherwise. ipiv is zero-indexed. // -// Getrf returns whether the matrix A is singular. The LU decomposition will -// be computed regardless of the singularity of A, but division by zero -// will occur if the false is returned and the result is used to solve a -// system of equations. +// Getrf returns whether the matrix A is nonsingular. The LU decomposition will +// be computed regardless of the singularity of A, but the result should not be +// used to solve a system of equation. func Getrf(a blas64.General, ipiv []int) bool { return lapack64.Dgetrf(a.Rows, a.Cols, a.Data, max(1, a.Stride), ipiv) } @@ -606,27 +642,50 @@ func Lantb(norm lapack.MatrixNorm, a blas64.TriangularBand, work []float64) floa // // X[i,0:n] is moved to X[k[i],0:n] for i=0,1,...,m-1. // -// k must have length m, otherwise Lapmr will panic. +// k must have length m, otherwise Lapmr will panic. k is zero-indexed. func Lapmr(forward bool, x blas64.General, k []int) { lapack64.Dlapmr(forward, x.Rows, x.Cols, x.Data, max(1, x.Stride), k) } // Lapmt rearranges the columns of the m×n matrix X as specified by the -// permutation k_0, k_1, ..., k_{n-1} of the integers 0, ..., n-1. +// permutation k[0],k[1],...,k[n-1] of the integers 0,...,n-1. // -// If forward is true a forward permutation is performed: +// If forward is true, a forward permutation is applied: // -// X[0:m, k[j]] is moved to X[0:m, j] for j = 0, 1, ..., n-1. +// X[0:m,k[j]] is moved to X[0:m,j] for j=0,1,...,n-1. // -// otherwise a backward permutation is performed: +// If forward is false, a backward permutation is applied: // -// X[0:m, j] is moved to X[0:m, k[j]] for j = 0, 1, ..., n-1. +// X[0:m,j] is moved to X[0:m,k[j]] for j=0,1,...,n-1. // // k must have length n, otherwise Lapmt will panic. k is zero-indexed. func Lapmt(forward bool, x blas64.General, k []int) { lapack64.Dlapmt(forward, x.Rows, x.Cols, x.Data, max(1, x.Stride), k) } +// Orglq generates an m×n matrix Q with orthonormal rows defined as the first m +// rows of a product of k elementary reflectors of order n +// +// Q = H_{k-1} * ... * H_0 +// +// as returned by Dgelqf. +// +// k is determined by the length of tau. +// +// On entry, tau and the first k rows of A must contain the scalar factors and +// the vectors, respectively, which define the elementary reflectors H_i, +// i=0,...,k-1, as returned by Dgelqf. On return, A contains the matrix Q. +// +// work must have length at least lwork and lwork must be at least max(1,m). On +// return, optimal value of lwork will be stored in work[0]. It must also hold +// that 0 <= k <= m <= n, otherwise Orglq will panic. +// +// If lwork == -1, instead of performing Orglq, the function only calculates the +// optimal value of lwork and stores it into work[0]. +func Orglq(a blas64.General, tau, work []float64, lwork int) { + lapack64.Dorglq(a.Rows, a.Cols, len(tau), a.Data, a.Stride, tau, work, lwork) +} + // Ormlq multiplies the matrix C by the othogonal matrix Q defined by // A and tau. A and tau are as returned from Gelqf. // @@ -651,6 +710,28 @@ func Ormlq(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64 lapack64.Dormlq(side, trans, c.Rows, c.Cols, a.Rows, a.Data, max(1, a.Stride), tau, c.Data, max(1, c.Stride), work, lwork) } +// Orgqr generates an m×n matrix Q with orthonormal columns defined by the +// product of elementary reflectors +// +// Q = H_0 * H_1 * ... * H_{k-1} +// +// as computed by Geqrf. +// +// k is determined by the length of tau. +// +// The length of work must be at least n and it also must be that 0 <= k <= n +// and 0 <= n <= m. +// +// work is temporary storage, and lwork specifies the usable memory length. At +// minimum, lwork >= n, and the amount of blocking is limited by the usable +// length. If lwork == -1, instead of computing Orgqr the optimal work length +// is stored into work[0]. +// +// Orgqr will panic if the conditions on input values are not met. +func Orgqr(a blas64.General, tau []float64, work []float64, lwork int) { + lapack64.Dorgqr(a.Rows, a.Cols, len(tau), a.Data, a.Stride, tau, work, lwork) +} + // Ormqr multiplies an m×n matrix C by an orthogonal matrix Q as // // C = Q * C if side == blas.Left and trans == blas.NoTrans, @@ -662,12 +743,13 @@ func Ormlq(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64 // // Q = H_0 * H_1 * ... * H_{k-1}. // +// k is determined by the length of tau. +// // If side == blas.Left, A is an m×k matrix and 0 <= k <= m. // If side == blas.Right, A is an n×k matrix and 0 <= k <= n. // The ith column of A contains the vector which defines the elementary -// reflector H_i and tau[i] contains its scalar factor. tau must have length k -// and Ormqr will panic otherwise. Geqrf returns A and tau in the required -// form. +// reflector H_i and tau[i] contains its scalar factor. Geqrf returns A and tau +// in the required form. // // work must have length at least max(1,lwork), and lwork must be at least n if // side == blas.Left and at least m if side == blas.Right, otherwise Ormqr will @@ -682,7 +764,7 @@ func Ormlq(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64 // If lwork is -1, instead of performing Ormqr, the optimal workspace size will // be stored into work[0]. func Ormqr(side blas.Side, trans blas.Transpose, a blas64.General, tau []float64, c blas64.General, work []float64, lwork int) { - lapack64.Dormqr(side, trans, c.Rows, c.Cols, a.Cols, a.Data, max(1, a.Stride), tau, c.Data, max(1, c.Stride), work, lwork) + lapack64.Dormqr(side, trans, c.Rows, c.Cols, len(tau), a.Data, max(1, a.Stride), tau, c.Data, max(1, c.Stride), work, lwork) } // Pocon estimates the reciprocal of the condition number of a positive-definite @@ -804,7 +886,7 @@ func Trtrs(trans blas.Transpose, a blas64.Triangular, b blas64.General) (ok bool // larger. On return, optimal value of lwork will be stored in work[0]. // // If lwork == -1, instead of performing Geev, the function only calculates the -// optimal vaule of lwork and stores it into work[0]. +// optimal value of lwork and stores it into work[0]. // // On return, first will be the index of the first valid eigenvalue. // If first == 0, all eigenvalues and eigenvectors have been computed. diff --git a/vendor/gonum.org/v1/gonum/mat/cholesky.go b/vendor/gonum.org/v1/gonum/mat/cholesky.go index 0f957cdd..f11948d0 100644 --- a/vendor/gonum.org/v1/gonum/mat/cholesky.go +++ b/vendor/gonum.org/v1/gonum/mat/cholesky.go @@ -25,6 +25,9 @@ var ( _ Symmetric = (*BandCholesky)(nil) _ Banded = (*BandCholesky)(nil) _ SymBanded = (*BandCholesky)(nil) + + _ Matrix = (*PivotedCholesky)(nil) + _ Symmetric = (*PivotedCholesky)(nil) ) // Cholesky is a symmetric positive definite matrix represented by its @@ -78,18 +81,12 @@ func (c *Cholesky) updateCond(norm float64) { // Dims returns the dimensions of the matrix. func (ch *Cholesky) Dims() (r, c int) { - if !ch.valid() { - panic(badCholesky) - } - r, c = ch.chol.Dims() - return r, c + n := ch.SymmetricDim() + return n, n } // At returns the element at row i, column j. func (c *Cholesky) At(i, j int) float64 { - if !c.valid() { - panic(badCholesky) - } n := c.SymmetricDim() if uint(i) >= uint(n) { panic(ErrRowAccess) @@ -113,8 +110,11 @@ func (c *Cholesky) T() Matrix { // SymmetricDim implements the Symmetric interface and returns the number of rows // in the matrix (this is also the number of columns). func (c *Cholesky) SymmetricDim() int { - r, _ := c.chol.Dims() - return r + if c.chol == nil { + return 0 + } + n, _ := c.chol.Triangle() + return n } // Cond returns the condition number of the factorized matrix. @@ -307,10 +307,15 @@ func (c *Cholesky) SolveVecTo(dst *VecDense, b Vector) error { } } -// RawU returns the Triangular matrix used to store the Cholesky decomposition of -// the original matrix A. The returned matrix should not be modified. If it is -// modified, the decomposition is invalid and should not be used. +// RawU returns the Triangular matrix used to store the Cholesky factorization +// of the original matrix A. If the returned matrix is modified, the +// factorization is invalid and should not be used. +// +// If Factorize has not been called, RawU will return nil. func (c *Cholesky) RawU() Triangular { + if !c.valid() { + return nil + } return c.chol } @@ -842,19 +847,13 @@ func (ch *BandCholesky) Reset() { // Dims returns the dimensions of the matrix. func (ch *BandCholesky) Dims() (r, c int) { - if !ch.valid() { - panic(badCholesky) - } - r, c = ch.chol.Dims() - return r, c + n := ch.SymmetricDim() + return n, n } // At returns the element at row i, column j. func (ch *BandCholesky) At(i, j int) float64 { - if !ch.valid() { - panic(badCholesky) - } - n, k, _ := ch.chol.TriBand() + n, k := ch.SymBand() if uint(i) >= uint(n) { panic(ErrRowAccess) } @@ -888,6 +887,9 @@ func (ch *BandCholesky) TBand() Banded { // SymmetricDim implements the Symmetric interface and returns the number of rows // in the matrix (this is also the number of columns). func (ch *BandCholesky) SymmetricDim() int { + if ch.chol == nil { + return 0 + } n, _ := ch.chol.Triangle() return n } @@ -936,3 +938,266 @@ func (ch *BandCholesky) LogDet() float64 { func (ch *BandCholesky) valid() bool { return ch.chol != nil && !ch.chol.IsEmpty() } + +// PivotedCholesky is a symmetric positive semi-definite matrix represented by +// its Cholesky factorization with complete pivoting. +// +// The factorization has the form +// +// Pᵀ * A * P = Uᵀ * U +// +// where U is an upper triangular matrix and P is a permutation matrix. +// +// Cholesky methods may only be called on a receiver that has been initialized +// by a call to Factorize. SolveTo and SolveVecTo methods may only called if +// Factorize has returned true. +// +// If the matrix A is certainly positive definite, then the unpivoted Cholesky +// could be more efficient, especially for smaller matrices. +type PivotedCholesky struct { + chol *TriDense // The factor U + piv, pivTrans []int // The permutation matrices P and Pᵀ + rank int // The computed rank of A + + ok bool // Indicates whether and the factorization can be used for solving linear systems + cond float64 // The condition number when ok is true +} + +// Factorize computes the Cholesky factorization of the symmetric positive +// semi-definite matrix A and returns whether the matrix is positive definite. +// If Factorize returns false, the SolveTo methods must not be used. +// +// tol is a tolerance used to determine the computed rank of A. If it is +// negative, a default value will be used. +func (c *PivotedCholesky) Factorize(a Symmetric, tol float64) (ok bool) { + n := a.SymmetricDim() + c.reset(n) + copySymIntoTriangle(c.chol, a) + + work := getFloat64s(3*c.chol.mat.N, false) + defer putFloat64s(work) + + sym := c.chol.asSymBlas() + aNorm := lapack64.Lansy(CondNorm, sym, work) + _, c.rank, c.ok = lapack64.Pstrf(sym, c.piv, tol, work) + if c.ok { + iwork := getInts(n, false) + defer putInts(iwork) + c.cond = 1 / lapack64.Pocon(sym, aNorm, work, iwork) + } else { + for i := c.rank; i < n; i++ { + zero(sym.Data[i*sym.Stride+i : i*sym.Stride+n]) + } + } + for i, p := range c.piv { + c.pivTrans[p] = i + } + + return c.ok +} + +// reset prepares the receiver for factorization of matrices of size n. +func (c *PivotedCholesky) reset(n int) { + if c.chol == nil { + c.chol = NewTriDense(n, Upper, nil) + } else { + c.chol.Reset() + c.chol.reuseAsNonZeroed(n, Upper) + } + c.piv = useInt(c.piv, n) + c.pivTrans = useInt(c.pivTrans, n) + c.rank = 0 + c.ok = false + c.cond = math.Inf(1) +} + +// Dims returns the dimensions of the matrix A. +func (ch *PivotedCholesky) Dims() (r, c int) { + n := ch.SymmetricDim() + return n, n +} + +// At returns the element of A at row i, column j. +func (c *PivotedCholesky) At(i, j int) float64 { + n := c.SymmetricDim() + if uint(i) >= uint(n) { + panic(ErrRowAccess) + } + if uint(j) >= uint(n) { + panic(ErrColAccess) + } + + i = c.pivTrans[i] + j = c.pivTrans[j] + minij := min(min(i+1, j+1), c.rank) + var val float64 + for k := 0; k < minij; k++ { + val += c.chol.at(k, i) * c.chol.at(k, j) + } + return val +} + +// T returns the receiver, the transpose of a symmetric matrix. +func (c *PivotedCholesky) T() Matrix { + return c +} + +// SymmetricDim implements the Symmetric interface and returns the number of +// rows (or columns) in the matrix . +func (c *PivotedCholesky) SymmetricDim() int { + if c.chol == nil { + return 0 + } + n, _ := c.chol.Triangle() + return n +} + +// Rank returns the computed rank of the matrix A. +func (c *PivotedCholesky) Rank() int { + if c.chol == nil { + panic(badCholesky) + } + return c.rank +} + +// Cond returns the condition number of the factorized matrix. +func (c *PivotedCholesky) Cond() float64 { + if c.chol == nil { + panic(badCholesky) + } + return c.cond +} + +// RawU returns the Triangular matrix used to store the Cholesky factorization +// of the original matrix A. If the returned matrix is modified, the +// factorization is invalid and should not be used. +// +// If Factorized returned false, the rows of U from Rank to n will contain zeros +// and so U will be upper trapezoidal. +// +// If Factorize has not been called, RawU will return nil. +func (c *PivotedCholesky) RawU() Triangular { + if c.chol == nil { + return nil + } + return c.chol +} + +// UTo stores the n×n upper triangular matrix U from the Cholesky factorization +// +// Pᵀ * A * P = Uᵀ * U. +// +// into dst. If dst is empty, it is resized to be an n×n upper triangular +// matrix. When dst is non-empty, UTo panics if dst is not n×n or not Upper. +// +// If Factorized returned false, the rows of U from Rank to n will contain zeros +// and so U will be upper trapezoidal. +func (c *PivotedCholesky) UTo(dst *TriDense) { + if c.chol == nil { + panic(badCholesky) + } + n := c.chol.mat.N + if dst.IsEmpty() { + dst.ReuseAsTri(n, Upper) + } else { + n2, kind := dst.Triangle() + if n != n2 { + panic(ErrShape) + } + if kind != Upper { + panic(ErrTriangle) + } + } + dst.Copy(c.chol) +} + +// ColumnPivots returns the column permutation p that represents the permutation +// matrix P from the Cholesky factorization +// +// Pᵀ * A * P = Uᵀ * U +// +// such that the nonzero entries are P[p[k],k] = 1. +func (c *PivotedCholesky) ColumnPivots(dst []int) []int { + if c.chol == nil { + panic(badCholesky) + } + n := c.chol.mat.N + if dst == nil { + dst = make([]int, n) + } + if len(dst) != n { + panic(badSliceLength) + } + copy(dst, c.piv) + return dst +} + +// SolveTo finds the matrix X that solves A * X = B where A is represented by +// the Cholesky decomposition. The result is stored in-place into dst. If the +// Cholesky decomposition is singular or near-singular, a Condition error is +// returned. See the documentation for Condition for more information. +// +// If Factorize returned false, SolveTo will panic. +func (c *PivotedCholesky) SolveTo(dst *Dense, b Matrix) error { + if !c.ok { + panic(badCholesky) + } + n := c.chol.mat.N + bm, bn := b.Dims() + if n != bm { + panic(ErrShape) + } + + dst.reuseAsNonZeroed(bm, bn) + if dst != b { + dst.Copy(b) + } + + // Permute rows of B: D = Pᵀ * B. + lapack64.Lapmr(true, dst.mat, c.piv) + // Solve Uᵀ * U * Y = D. + lapack64.Potrs(c.chol.mat, dst.mat) + // Permute rows of Y to recover the solution: X = P * Y. + lapack64.Lapmr(false, dst.mat, c.piv) + + if c.cond > ConditionTolerance { + return Condition(c.cond) + } + return nil +} + +// SolveVecTo finds the vector x that solves A * x = b where A is represented by +// the Cholesky decomposition. The result is stored in-place into dst. If the +// Cholesky decomposition is singular or near-singular, a Condition error is +// returned. See the documentation for Condition for more information. +// +// If Factorize returned false, SolveVecTo will panic. +func (c *PivotedCholesky) SolveVecTo(dst *VecDense, b Vector) error { + if !c.ok { + panic(badCholesky) + } + n := c.chol.mat.N + if br, bc := b.Dims(); br != n || bc != 1 { + panic(ErrShape) + } + if b, ok := b.(RawVectorer); ok && dst != b { + dst.checkOverlap(b.RawVector()) + } + + dst.reuseAsNonZeroed(n) + if dst != b { + dst.CopyVec(b) + } + + // Permute rows of B: D = Pᵀ * B. + lapack64.Lapmr(true, dst.asGeneral(), c.piv) + // Solve Uᵀ * U * Y = D. + lapack64.Potrs(c.chol.mat, dst.asGeneral()) + // Permute rows of Y to recover the solution: X = P * Y. + lapack64.Lapmr(false, dst.asGeneral(), c.piv) + + if c.cond > ConditionTolerance { + return Condition(c.cond) + } + return nil +} diff --git a/vendor/gonum.org/v1/gonum/mat/dense.go b/vendor/gonum.org/v1/gonum/mat/dense.go index 5ec425fa..b08360cc 100644 --- a/vendor/gonum.org/v1/gonum/mat/dense.go +++ b/vendor/gonum.org/v1/gonum/mat/dense.go @@ -613,3 +613,58 @@ func (m *Dense) Norm(norm float64) float64 { } return lapack64.Lange(lnorm, m.mat, nil) } + +// Permutation constructs an n×n permutation matrix P from the given +// row permutation such that the nonzero entries are P[i,p[i]] = 1. +func (m *Dense) Permutation(n int, p []int) { + if len(p) != n { + panic(badSliceLength) + } + m.reuseAsZeroed(n, n) + for i, v := range p { + if v < 0 || v >= n { + panic(ErrRowAccess) + } + m.mat.Data[i*m.mat.Stride+v] = 1 + } +} + +// PermuteRows rearranges the rows of the m×n matrix A in the receiver as +// specified by the permutation p[0],p[1],...,p[m-1] of the integers 0,...,m-1. +// +// If inverse is false, the given permutation is applied: +// +// A[p[i],0:n] is moved to A[i,0:n] for i=0,1,...,m-1. +// +// If inverse is true, the inverse permutation is applied: +// +// A[i,0:n] is moved to A[p[i],0:n] for i=0,1,...,m-1. +// +// p must have length m, otherwise PermuteRows will panic. +func (m *Dense) PermuteRows(p []int, inverse bool) { + r, _ := m.Dims() + if len(p) != r { + panic(badSliceLength) + } + lapack64.Lapmr(!inverse, m.mat, p) +} + +// PermuteCols rearranges the columns of the m×n matrix A in the reciever as +// specified by the permutation p[0],p[1],...,p[n-1] of the integers 0,...,n-1. +// +// If inverse is false, the given permutation is applied: +// +// A[0:m,p[j]] is moved to A[0:m,j] for j = 0, 1, ..., n-1. +// +// If inverse is true, the inverse permutation is applied: +// +// A[0:m,j] is moved to A[0:m,p[j]] for j = 0, 1, ..., n-1. +// +// p must have length n, otherwise PermuteCols will panic. +func (m *Dense) PermuteCols(p []int, inverse bool) { + _, c := m.Dims() + if len(p) != c { + panic(badSliceLength) + } + lapack64.Lapmt(!inverse, m.mat, p) +} diff --git a/vendor/gonum.org/v1/gonum/mat/eigen.go b/vendor/gonum.org/v1/gonum/mat/eigen.go index 69f8eb57..859247d8 100644 --- a/vendor/gonum.org/v1/gonum/mat/eigen.go +++ b/vendor/gonum.org/v1/gonum/mat/eigen.go @@ -14,8 +14,12 @@ const ( noVectors = "mat: eigenvectors not computed" ) -// EigenSym is a type for creating and manipulating the Eigen decomposition of -// symmetric matrices. +// EigenSym is a type for computing all eigenvalues and, optionally, +// eigenvectors of a symmetric matrix A. +// +// It is a Symmetric matrix represented by its spectral factorization. Once +// computed, this representation is useful for extracting eigenvalues and +// eigenvector, but At is slow. type EigenSym struct { vectorsComputed bool @@ -23,18 +27,59 @@ type EigenSym struct { vectors *Dense } -// Factorize computes the eigenvalue decomposition of the symmetric matrix a. -// The Eigen decomposition is defined as +// Dims returns the dimensions of the matrix. +func (e *EigenSym) Dims() (r, c int) { + n := e.SymmetricDim() + return n, n +} + +// SymmetricDim implements the Symmetric interface. +func (e *EigenSym) SymmetricDim() int { + return len(e.values) +} + +// At returns the element at row i, column j of the matrix A. +// +// At will panic if the eigenvectors have not been computed. +func (e *EigenSym) At(i, j int) float64 { + if !e.vectorsComputed { + panic(noVectors) + } + n, _ := e.Dims() + if uint(i) >= uint(n) { + panic(ErrRowAccess) + } + if uint(j) >= uint(n) { + panic(ErrColAccess) + } + + var val float64 + for k := 0; k < n; k++ { + val += e.values[k] * e.vectors.at(i, k) * e.vectors.at(j, k) + } + return val +} + +// T returns the receiver, the transpose of a symmetric matrix. +func (e *EigenSym) T() Matrix { + return e +} + +// Factorize computes the spectral factorization (eigendecomposition) of the +// symmetric matrix A. // -// A = P * D * P^-1 +// The spectral factorization of A can be written as // -// where D is a diagonal matrix containing the eigenvalues of the matrix, and -// P is a matrix of the eigenvectors of A. Factorize computes the eigenvalues -// in ascending order. If the vectors input argument is false, the eigenvectors -// are not computed. +// A = Q * Λ * Qᵀ // -// Factorize returns whether the decomposition succeeded. If the decomposition -// failed, methods that require a successful factorization will panic. +// where Λ is a diagonal matrix whose entries are the eigenvalues, and Q is an +// orthogonal matrix whose columns are the eigenvectors. +// +// If vectors is false, the eigenvectors are not computed and later calls to +// VectorsTo and At will panic. +// +// Factorize returns whether the factorization succeeded. If it returns false, +// methods that require a successful factorization will panic. func (e *EigenSym) Factorize(a Symmetric, vectors bool) (ok bool) { // kill previous decomposition e.vectorsComputed = false @@ -72,12 +117,15 @@ func (e *EigenSym) succFact() bool { return len(e.values) != 0 } -// Values extracts the eigenvalues of the factorized matrix in ascending order. -// If dst is 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 nil, then a new -// slice will be allocated of the proper length and filled with the eigenvalues. +// Values extracts the eigenvalues of the factorized n×n matrix A in ascending +// order. // -// Values panics if the Eigen decomposition was not successful. +// If dst is not nil, the values are stored in-place into dst and returned, +// otherwise a new slice is allocated first. If dst is not nil, it must have +// length equal to n. +// +// If the receiver does not contain a successful factorization, Values will +// panic. func (e *EigenSym) Values(dst []float64) []float64 { if !e.succFact() { panic(badFact) @@ -92,13 +140,27 @@ func (e *EigenSym) Values(dst []float64) []float64 { return dst } -// VectorsTo stores the eigenvectors of the decomposition into the columns of -// dst. +// RawValues returns the slice storing the eigenvalues of A in ascending order. // -// If dst is empty, VectorsTo will resize dst to be n×n. When dst is -// non-empty, VectorsTo will panic if dst is not n×n. VectorsTo will also -// panic if the eigenvectors were not computed during the factorization, -// or if the receiver does not contain a successful factorization. +// If the returned slice is modified, the factorization is invalid and should +// not be used. +// +// If the receiver does not contain a successful factorization, RawValues will +// return nil. +func (e *EigenSym) RawValues() []float64 { + if !e.succFact() { + return nil + } + return e.values +} + +// VectorsTo stores the orthonormal eigenvectors of the factorized n×n matrix A +// into the columns of dst. +// +// If dst is empty, VectorsTo will resize dst to be n×n. When dst is non-empty, +// VectorsTo will panic if dst is not n×n. VectorsTo will also panic if the +// eigenvectors were not computed during the factorization, or if the receiver +// does not contain a successful factorization. func (e *EigenSym) VectorsTo(dst *Dense) { if !e.succFact() { panic(badFact) @@ -118,6 +180,25 @@ func (e *EigenSym) VectorsTo(dst *Dense) { dst.Copy(e.vectors) } +// RawQ returns the orthogonal matrix Q from the spectral factorization of the +// original matrix A +// +// A = Q * Λ * Qᵀ +// +// The columns of Q contain the eigenvectors of A. +// +// If the returned matrix is modified, the factorization is invalid and should +// not be used. +// +// If the receiver does not contain a successful factorization or eigenvectors +// not computed, RawU will return nil. +func (e *EigenSym) RawQ() Matrix { + if !e.succFact() || !e.vectorsComputed { + return nil + } + return e.vectors +} + // EigenKind specifies the computation of eigenvectors during factorization. type EigenKind int diff --git a/vendor/gonum.org/v1/gonum/mat/lq.go b/vendor/gonum.org/v1/gonum/mat/lq.go index 4d6da5ee..a3b3543b 100644 --- a/vendor/gonum.org/v1/gonum/mat/lq.go +++ b/vendor/gonum.org/v1/gonum/mat/lq.go @@ -18,10 +18,42 @@ const badLQ = "mat: invalid LQ factorization" // LQ is a type for creating and using the LQ factorization of a matrix. type LQ struct { lq *Dense + q *Dense tau []float64 cond float64 } +// Dims returns the dimensions of the matrix. +func (lq *LQ) Dims() (r, c int) { + if lq.lq == nil { + return 0, 0 + } + return lq.lq.Dims() +} + +// At returns the element at row i, column j. +func (lq *LQ) At(i, j int) float64 { + m, n := lq.Dims() + if uint(i) >= uint(m) { + panic(ErrRowAccess) + } + if uint(j) >= uint(n) { + panic(ErrColAccess) + } + + var val float64 + for k := 0; k <= i; k++ { + val += lq.lq.at(i, k) * lq.q.at(k, j) + } + return val +} + +// T performs an implicit transpose by returning the receiver inside a +// Transpose. +func (lq *LQ) T() Matrix { + return Transpose{lq} +} + func (lq *LQ) updateCond(norm lapack.MatrixNorm) { // Since A = L*Q, and Q is orthogonal, we get for the condition number κ // κ(A) := |A| |A^-1| = |L*Q| |(L*Q)^-1| = |L| |Qᵀ * L^-1| @@ -55,18 +87,34 @@ func (lq *LQ) factorize(a Matrix, norm lapack.MatrixNorm) { if m > n { panic(ErrShape) } - k := min(m, n) if lq.lq == nil { lq.lq = &Dense{} } lq.lq.CloneFrom(a) work := []float64{0} - lq.tau = make([]float64, k) + lq.tau = make([]float64, m) lapack64.Gelqf(lq.lq.mat, lq.tau, work, -1) work = getFloat64s(int(work[0]), false) lapack64.Gelqf(lq.lq.mat, lq.tau, work, len(work)) putFloat64s(work) lq.updateCond(norm) + lq.updateQ() +} + +func (lq *LQ) updateQ() { + _, n := lq.Dims() + if lq.q == nil { + lq.q = NewDense(n, n, nil) + } else { + lq.q.reuseAsNonZeroed(n, n) + } + // Construct Q from the elementary reflectors. + lq.q.Copy(lq.lq) + work := []float64{0} + lapack64.Orglq(lq.q.mat, lq.tau, work, -1) + work = getFloat64s(int(work[0]), false) + lapack64.Orglq(lq.q.mat, lq.tau, work, len(work)) + putFloat64s(work) } // isValid returns whether the receiver contains a factorization. @@ -130,38 +178,24 @@ func (lq *LQ) LTo(dst *Dense) { // QTo extracts the n×n orthonormal matrix Q from an LQ decomposition. // -// If dst is empty, QTo will resize dst to be c×c. When dst is -// non-empty, QTo will panic if dst is not c×c. QTo will also panic +// If dst is empty, QTo will resize dst to be n×n. When dst is +// non-empty, QTo will panic if dst is not n×n. QTo will also panic // if the receiver does not contain a successful factorization. func (lq *LQ) QTo(dst *Dense) { if !lq.isValid() { panic(badLQ) } - _, c := lq.lq.Dims() + _, n := lq.lq.Dims() if dst.IsEmpty() { - dst.ReuseAs(c, c) + dst.ReuseAs(n, n) } else { - r2, c2 := dst.Dims() - if c != r2 || c != c2 { + m2, n2 := dst.Dims() + if n != m2 || n != n2 { panic(ErrShape) } - dst.Zero() } - q := dst.mat - - // Set Q = I. - ldq := q.Stride - for i := 0; i < c; i++ { - q.Data[i*ldq+i] = 1 - } - - // Construct Q from the elementary reflectors. - work := []float64{0} - lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, -1) - work = getFloat64s(int(work[0]), false) - lapack64.Ormlq(blas.Left, blas.NoTrans, lq.lq.mat, lq.tau, q, work, len(work)) - putFloat64s(work) + dst.Copy(lq.q) } // SolveTo finds a minimum-norm solution to a system of linear equations defined diff --git a/vendor/gonum.org/v1/gonum/mat/lu.go b/vendor/gonum.org/v1/gonum/mat/lu.go index 4ec424f9..b530ada7 100644 --- a/vendor/gonum.org/v1/gonum/mat/lu.go +++ b/vendor/gonum.org/v1/gonum/mat/lu.go @@ -19,11 +19,62 @@ const ( badLU = "mat: invalid LU factorization" ) -// LU is a type for creating and using the LU factorization of a matrix. +// LU is a square n×n matrix represented by its LU factorization with partial +// pivoting. +// +// The factorization has the form +// +// A = P * L * U +// +// where P is a permutation matrix, L is lower triangular with unit diagonal +// elements, and U is upper triangular. +// +// Note that this matrix representation is useful for certain operations, in +// particular for solving linear systems of equations. It is very inefficient at +// other operations, in particular At is slow. type LU struct { lu *Dense - pivot []int + swaps []int + piv []int cond float64 + ok bool // Whether A is nonsingular +} + +var _ Matrix = (*LU)(nil) + +// Dims returns the dimensions of the matrix A. +func (lu *LU) Dims() (r, c int) { + if lu.lu == nil { + return 0, 0 + } + return lu.lu.Dims() +} + +// At returns the element of A at row i, column j. +func (lu *LU) At(i, j int) float64 { + n, _ := lu.Dims() + if uint(i) >= uint(n) { + panic(ErrRowAccess) + } + if uint(j) >= uint(n) { + panic(ErrColAccess) + } + + i = lu.piv[i] + var val float64 + for k := 0; k < min(i, j+1); k++ { + val += lu.lu.at(i, k) * lu.lu.at(k, j) + } + if i <= j { + val += lu.lu.at(i, j) + } + return val +} + +// T performs an implicit transpose by returning the receiver inside a +// Transpose. +func (lu *LU) T() Matrix { + return Transpose{lu} } // updateCond updates the stored condition number of the matrix. anorm is the @@ -52,40 +103,51 @@ func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) { lu.cond = 1 / v } -// Factorize computes the LU factorization of the square matrix a and stores the -// result. The LU decomposition will complete regardless of the singularity of a. +// Factorize computes the LU factorization of the square matrix A and stores the +// result in the receiver. The LU decomposition will complete regardless of the +// singularity of a. // -// The LU factorization is computed with pivoting, and so really the decomposition -// is a PLU decomposition where P is a permutation matrix. The individual matrix -// factors can be extracted from the factorization using the Permutation method -// on Dense, and the LU.LTo and LU.UTo methods. +// The L and U matrix factors can be extracted from the factorization using the +// LTo and UTo methods. The matrix P can be extracted as a row permutation using +// the RowPivots method and applied using Dense.PermuteRows. func (lu *LU) Factorize(a Matrix) { lu.factorize(a, CondNorm) } func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) { - r, c := a.Dims() - if r != c { + m, n := a.Dims() + if m != n { panic(ErrSquare) } if lu.lu == nil { - lu.lu = NewDense(r, r, nil) + lu.lu = NewDense(n, n, nil) } else { lu.lu.Reset() - lu.lu.reuseAsNonZeroed(r, r) + lu.lu.reuseAsNonZeroed(n, n) } lu.lu.Copy(a) - if cap(lu.pivot) < r { - lu.pivot = make([]int, r) - } - lu.pivot = lu.pivot[:r] - work := getFloat64s(r, false) + lu.swaps = useInt(lu.swaps, n) + lu.piv = useInt(lu.piv, n) + work := getFloat64s(n, false) anorm := lapack64.Lange(norm, lu.lu.mat, work) putFloat64s(work) - lapack64.Getrf(lu.lu.mat, lu.pivot) + lu.ok = lapack64.Getrf(lu.lu.mat, lu.swaps) + lu.updatePivots(lu.swaps) lu.updateCond(anorm, norm) } +func (lu *LU) updatePivots(swaps []int) { + // Replay the sequence of row swaps in order to find the row permutation. + for i := range lu.piv { + lu.piv[i] = i + } + n, _ := lu.Dims() + for i := n - 1; i >= 0; i-- { + v := swaps[i] + lu.piv[i], lu.piv[v] = lu.piv[v], lu.piv[i] + } +} + // isValid returns whether the receiver contains a factorization. func (lu *LU) isValid() bool { return lu.lu != nil && !lu.lu.IsEmpty() @@ -106,17 +168,21 @@ func (lu *LU) Reset() { if lu.lu != nil { lu.lu.Reset() } - lu.pivot = lu.pivot[:0] + lu.swaps = lu.swaps[:0] + lu.piv = lu.piv[:0] } func (lu *LU) isZero() bool { - return len(lu.pivot) == 0 + return len(lu.swaps) == 0 } // Det returns the determinant of the matrix that has been factorized. In many // expressions, using LogDet will be more numerically stable. // Det will panic if the receiver does not contain a factorization. func (lu *LU) Det() float64 { + if !lu.ok { + return 0 + } det, sign := lu.LogDet() return math.Exp(det) * sign } @@ -139,7 +205,7 @@ func (lu *LU) LogDet() (det float64, sign float64) { if v < 0 { sign *= -1 } - if lu.pivot[i] != i { + if lu.swaps[i] != i { sign *= -1 } logDiag[i] = math.Log(math.Abs(v)) @@ -147,39 +213,41 @@ func (lu *LU) LogDet() (det float64, sign float64) { return floats.Sum(logDiag), sign } -// Pivot returns pivot indices that enable the construction of the permutation -// matrix P (see Dense.Permutation). If swaps == nil, then new memory will be -// allocated, otherwise the length of the input must be equal to the size of the -// factorized matrix. -// Pivot will panic if the receiver does not contain a factorization. -func (lu *LU) Pivot(swaps []int) []int { +// RowPivots returns the row permutation that represents the permutation matrix +// P from the LU factorization +// +// A = P * L * U. +// +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal the size of the factorized matrix, RowPivots +// will panic. RowPivots will panic if the receiver does not contain a +// factorization. +func (lu *LU) RowPivots(dst []int) []int { if !lu.isValid() { panic(badLU) } - _, n := lu.lu.Dims() - if swaps == nil { - swaps = make([]int, n) + if dst == nil { + dst = make([]int, n) } - if len(swaps) != n { + if len(dst) != n { panic(badSliceLength) } - // Perform the inverse of the row swaps in order to find the final - // row swap position. - for i := range swaps { - swaps[i] = i - } - for i := n - 1; i >= 0; i-- { - v := lu.pivot[i] - swaps[i], swaps[v] = swaps[v], swaps[i] - } - return swaps + copy(dst, lu.piv) + return dst +} + +// Pivot returns the row pivots of the receiver. +// +// Deprecated: Use RowPivots instead. +func (lu *LU) Pivot(dst []int) []int { + return lu.RowPivots(dst) } // RankOne updates an LU factorization as if a rank-one update had been applied to // the original matrix A, storing the result into the receiver. That is, if in // the original LU decomposition P * L * U = A, in the updated decomposition -// P * L * U = A + alpha * x * yᵀ. +// P * L' * U' = A + alpha * x * yᵀ. // RankOne will panic if orig does not contain a factorization. func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) { if !orig.isValid() { @@ -198,19 +266,18 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) { } if orig != lu { if lu.isZero() { - if cap(lu.pivot) < n { - lu.pivot = make([]int, n) - } - lu.pivot = lu.pivot[:n] + lu.swaps = useInt(lu.swaps, n) + lu.piv = useInt(lu.piv, n) if lu.lu == nil { lu.lu = NewDense(n, n, nil) } else { lu.lu.reuseAsNonZeroed(n, n) } - } else if len(lu.pivot) != n { + } else if len(lu.swaps) != n { panic(ErrShape) } - copy(lu.pivot, orig.pivot) + copy(lu.swaps, orig.swaps) + lu.updatePivots(lu.swaps) lu.lu.Copy(orig.lu) } @@ -224,7 +291,7 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) { } // Adjust for the pivoting in the LU factorization - for i, v := range lu.pivot { + for i, v := range lu.swaps { xs[i], xs[v] = xs[v], xs[i] } @@ -273,10 +340,8 @@ func (lu *LU) LTo(dst *TriDense) *TriDense { } } // Extract the lower triangular elements. - for i := 0; i < n; i++ { - for j := 0; j < i; j++ { - dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] - } + for i := 1; i < n; i++ { + copy(dst.mat.Data[i*dst.mat.Stride:i*dst.mat.Stride+i], lu.lu.mat.Data[i*lu.lu.mat.Stride:i*lu.lu.mat.Stride+i]) } // Set ones on the diagonal. for i := 0; i < n; i++ { @@ -310,40 +375,21 @@ func (lu *LU) UTo(dst *TriDense) { } // Extract the upper triangular elements. for i := 0; i < n; i++ { - for j := i; j < n; j++ { - dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] - } - } -} - -// Permutation constructs an r×r permutation matrix with the given row swaps. -// A permutation matrix has exactly one element equal to one in each row and column -// and all other elements equal to zero. swaps[i] specifies the row with which -// i will be swapped, which is equivalent to the non-zero column of row i. -func (m *Dense) Permutation(r int, swaps []int) { - m.reuseAsNonZeroed(r, r) - for i := 0; i < r; i++ { - zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+r]) - v := swaps[i] - if v < 0 || v >= r { - panic(ErrRowAccess) - } - m.mat.Data[i*m.mat.Stride+v] = 1 + copy(dst.mat.Data[i*dst.mat.Stride+i:i*dst.mat.Stride+n], lu.lu.mat.Data[i*lu.lu.mat.Stride+i:i*lu.lu.mat.Stride+n]) } } -// SolveTo solves a system of linear equations using the LU decomposition of a matrix. -// It computes +// SolveTo solves a system of linear equations // -// A * X = B if trans == false -// Aᵀ * X = B if trans == true +// A * X = B if trans == false +// Aᵀ * X = B if trans == true // -// In both cases, A is represented in LU factorized form, and the matrix X is -// stored into dst. +// using the LU factorization of A stored in the receiver. The solution matrix X +// is stored into dst. // -// If A is singular or near-singular a Condition error is returned. See -// the documentation for Condition for more information. -// SolveTo will panic if the receiver does not contain a factorization. +// If A is singular or near-singular a Condition error is returned. See the +// documentation for Condition for more information. SolveTo will panic if the +// receiver does not contain a factorization. func (lu *LU) SolveTo(dst *Dense, trans bool, b Matrix) error { if !lu.isValid() { panic(badLU) @@ -354,16 +400,15 @@ func (lu *LU) SolveTo(dst *Dense, trans bool, b Matrix) error { if br != n { panic(ErrShape) } - // TODO(btracey): Should test the condition number instead of testing that - // the determinant is exactly zero. - if lu.Det() == 0 { + + if !lu.ok { return Condition(math.Inf(1)) } dst.reuseAsNonZeroed(n, bc) bU, _ := untranspose(b) - var restore func() if dst == bU { + var restore func() dst, restore = dst.isolatedWorkspace(bU) defer restore() } else if rm, ok := bU.(RawMatrixer); ok { @@ -375,25 +420,24 @@ func (lu *LU) SolveTo(dst *Dense, trans bool, b Matrix) error { if trans { t = blas.Trans } - lapack64.Getrs(t, lu.lu.mat, dst.mat, lu.pivot) + lapack64.Getrs(t, lu.lu.mat, dst.mat, lu.swaps) if lu.cond > ConditionTolerance { return Condition(lu.cond) } return nil } -// SolveVecTo solves a system of linear equations using the LU decomposition of a matrix. -// It computes +// SolveVecTo solves a system of linear equations // -// A * x = b if trans == false -// Aᵀ * x = b if trans == true +// A * x = b if trans == false +// Aᵀ * x = b if trans == true // -// In both cases, A is represented in LU factorized form, and the vector x is -// stored into dst. +// using the LU factorization of A stored in the receiver. The solution matrix x +// is stored into dst. // -// If A is singular or near-singular a Condition error is returned. See -// the documentation for Condition for more information. -// SolveVecTo will panic if the receiver does not contain a factorization. +// If A is singular or near-singular a Condition error is returned. See the +// documentation for Condition for more information. SolveVecTo will panic if the +// receiver does not contain a factorization. func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error { if !lu.isValid() { panic(badLU) @@ -403,6 +447,7 @@ func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error { if br, bc := b.Dims(); br != n || bc != 1 { panic(ErrShape) } + switch rv := b.(type) { default: dst.reuseAsNonZeroed(n) @@ -411,9 +456,8 @@ func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error { if dst != b { dst.checkOverlap(rv.RawVector()) } - // TODO(btracey): Should test the condition number instead of testing that - // the determinant is exactly zero. - if lu.Det() == 0 { + + if !lu.ok { return Condition(math.Inf(1)) } @@ -434,7 +478,7 @@ func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error { if trans { t = blas.Trans } - lapack64.Getrs(t, lu.lu.mat, vMat, lu.pivot) + lapack64.Getrs(t, lu.lu.mat, vMat, lu.swaps) if lu.cond > ConditionTolerance { return Condition(lu.cond) } diff --git a/vendor/gonum.org/v1/gonum/mat/matrix.go b/vendor/gonum.org/v1/gonum/mat/matrix.go index 51093aaa..2d67bbe0 100644 --- a/vendor/gonum.org/v1/gonum/mat/matrix.go +++ b/vendor/gonum.org/v1/gonum/mat/matrix.go @@ -219,6 +219,17 @@ type ColNonZeroDoer interface { DoColNonZero(j int, fn func(i, j int, v float64)) } +// A SolveToer can solve a linear system A⋅X = B or Aᵀ⋅X = B where A is a matrix +// represented by the receiver and B is a given matrix, storing the result into +// dst. +// +// If dst is empty, SolveTo will resize it to the correct size, otherwise it +// must have the correct size. Individual implementations may impose other +// restrictions on the input parameters, for example that A is a square matrix. +type SolveToer interface { + SolveTo(dst *Dense, trans bool, b Matrix) error +} + // untranspose untransposes a matrix if applicable. If a is an Untransposer, then // untranspose returns the underlying matrix and true. If it is not, then it returns // the input matrix and false. @@ -951,20 +962,6 @@ func Trace(a Matrix) float64 { return v } -func min(a, b int) int { - if a < b { - return a - } - return b -} - -func max(a, b int) int { - if a > b { - return a - } - return b -} - // use returns a float64 slice with l elements, using f if it // has the necessary capacity, otherwise creating a new slice. func use(f []float64, l int) []float64 { diff --git a/vendor/gonum.org/v1/gonum/mat/qr.go b/vendor/gonum.org/v1/gonum/mat/qr.go index f54bcc86..7f8fec8f 100644 --- a/vendor/gonum.org/v1/gonum/mat/qr.go +++ b/vendor/gonum.org/v1/gonum/mat/qr.go @@ -18,10 +18,80 @@ const badQR = "mat: invalid QR factorization" // QR is a type for creating and using the QR factorization of a matrix. type QR struct { qr *Dense + q *Dense tau []float64 cond float64 } +// Dims returns the dimensions of the matrix. +func (qr *QR) Dims() (r, c int) { + if qr.qr == nil { + return 0, 0 + } + return qr.qr.Dims() +} + +// At returns the element at row i, column j. At will panic if the receiver +// does not contain a successful factorization. +func (qr *QR) At(i, j int) float64 { + if !qr.isValid() { + panic(badQR) + } + + m, n := qr.Dims() + if uint(i) >= uint(m) { + panic(ErrRowAccess) + } + if uint(j) >= uint(n) { + panic(ErrColAccess) + } + + if qr.q == nil || qr.q.IsEmpty() { + // Calculate Qi, Q i-th row + qi := getFloat64s(m, true) + qr.qRowTo(i, qi) + + // Compute QR(i,j) + var val float64 + for k := 0; k <= j; k++ { + val += qi[k] * qr.qr.at(k, j) + } + putFloat64s(qi) + return val + } + + var val float64 + for k := 0; k <= j; k++ { + val += qr.q.at(i, k) * qr.qr.at(k, j) + } + return val +} + +// qRowTo extracts the i-th row of the orthonormal matrix Q from a QR +// decomposition. +func (qr *QR) qRowTo(i int, dst []float64) { + c := blas64.General{ + Rows: 1, + Cols: len(dst), + Stride: len(dst), + Data: dst, + } + c.Data[i] = 1 // C is the i-th unit vector + + // Construct Qi from the elementary reflectors: Qi = C * (H(1) H(2) ... H(nTau)) + work := []float64{0} + lapack64.Ormqr(blas.Right, blas.NoTrans, qr.qr.mat, qr.tau, c, work, -1) + work = getFloat64s(int(work[0]), false) + lapack64.Ormqr(blas.Right, blas.NoTrans, qr.qr.mat, qr.tau, c, work, len(work)) + putFloat64s(work) +} + +// T performs an implicit transpose by returning the receiver inside a +// Transpose. +func (qr *QR) T() Matrix { + return Transpose{qr} +} + func (qr *QR) updateCond(norm lapack.MatrixNorm) { // Since A = Q*R, and Q is orthogonal, we get for the condition number κ // κ(A) := |A| |A^-1| = |Q*R| |(Q*R)^-1| = |R| |R^-1 * Qᵀ| @@ -55,18 +125,36 @@ func (qr *QR) factorize(a Matrix, norm lapack.MatrixNorm) { if m < n { panic(ErrShape) } - k := min(m, n) if qr.qr == nil { qr.qr = &Dense{} } qr.qr.CloneFrom(a) work := []float64{0} - qr.tau = make([]float64, k) + qr.tau = make([]float64, n) lapack64.Geqrf(qr.qr.mat, qr.tau, work, -1) work = getFloat64s(int(work[0]), false) lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work)) putFloat64s(work) qr.updateCond(norm) + if qr.q != nil { + qr.q.Reset() + } +} + +func (qr *QR) updateQ() { + m, _ := qr.Dims() + if qr.q == nil { + qr.q = NewDense(m, m, nil) + } else { + qr.q.reuseAsNonZeroed(m, m) + } + // Construct Q from the elementary reflectors. + qr.q.Copy(qr.qr) + work := []float64{0} + lapack64.Orgqr(qr.q.mat, qr.tau, work, -1) + work = getFloat64s(int(work[0]), false) + lapack64.Orgqr(qr.q.mat, qr.tau, work, len(work)) + putFloat64s(work) } // isValid returns whether the receiver contains a factorization. @@ -101,7 +189,7 @@ func (qr *QR) RTo(dst *Dense) { dst.ReuseAs(r, c) } else { r2, c2 := dst.Dims() - if c != r2 || c != c2 { + if r != r2 || c != c2 { panic(ErrShape) } } @@ -143,20 +231,12 @@ func (qr *QR) QTo(dst *Dense) { if r != r2 || r != c2 { panic(ErrShape) } - dst.Zero() } - // Set Q = I. - for i := 0; i < r*r; i += r + 1 { - dst.mat.Data[i] = 1 + if qr.q == nil || qr.q.IsEmpty() { + qr.updateQ() } - - // Construct Q from the elementary reflectors. - work := []float64{0} - lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1) - work = getFloat64s(int(work[0]), false) - lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work)) - putFloat64s(work) + dst.Copy(qr.q) } // SolveTo finds a minimum-norm solution to a system of linear equations defined diff --git a/vendor/gonum.org/v1/gonum/mat/solve.go b/vendor/gonum.org/v1/gonum/mat/solve.go index 7b7f3f6b..ffccce8c 100644 --- a/vendor/gonum.org/v1/gonum/mat/solve.go +++ b/vendor/gonum.org/v1/gonum/mat/solve.go @@ -4,12 +4,6 @@ package mat -import ( - "gonum.org/v1/gonum/blas" - "gonum.org/v1/gonum/blas/blas64" - "gonum.org/v1/gonum/lapack/lapack64" -) - // Solve solves the linear least squares problem // // minimize over x |b - A*x|_2 @@ -30,9 +24,17 @@ import ( // single call. Vectors b are stored in the columns of the m×k matrix B. Vectors // x will be stored in-place into the n×k receiver. // +// If the underlying matrix of a is a SolveToer, its SolveTo method is used, +// otherwise a Dense copy of a will be used for the solution. +// // If A does not have full rank, a Condition error is returned. See the // documentation for Condition for more information. func (m *Dense) Solve(a, b Matrix) error { + aU, aTrans := untransposeExtract(a) + if a, ok := aU.(SolveToer); ok { + return a.SolveTo(m, aTrans, b) + } + ar, ac := a.Dims() br, bc := b.Dims() if ar != br { @@ -40,54 +42,6 @@ func (m *Dense) Solve(a, b Matrix) error { } m.reuseAsNonZeroed(ac, bc) - // TODO(btracey): Add special cases for SymDense, etc. - aU, aTrans := untranspose(a) - bU, bTrans := untranspose(b) - switch rma := aU.(type) { - case RawTriangular: - side := blas.Left - tA := blas.NoTrans - if aTrans { - tA = blas.Trans - } - - switch rm := bU.(type) { - case RawMatrixer: - if m != bU || bTrans { - if m == bU || m.checkOverlap(rm.RawMatrix()) { - tmp := getDenseWorkspace(br, bc, false) - tmp.Copy(b) - m.Copy(tmp) - putDenseWorkspace(tmp) - break - } - m.Copy(b) - } - default: - if m != bU { - m.Copy(b) - } else if bTrans { - // m and b share data so Copy cannot be used directly. - tmp := getDenseWorkspace(br, bc, false) - tmp.Copy(b) - m.Copy(tmp) - putDenseWorkspace(tmp) - } - } - - rm := rma.RawTriangular() - blas64.Trsm(side, tA, 1, rm, m.mat) - work := getFloat64s(3*rm.N, false) - iwork := getInts(rm.N, false) - cond := lapack64.Trcon(CondNorm, rm, work, iwork) - putFloat64s(work) - putInts(iwork) - if cond > ConditionTolerance { - return Condition(cond) - } - return nil - } - switch { case ar == ac: if a == b { diff --git a/vendor/gonum.org/v1/gonum/mat/triangular.go b/vendor/gonum.org/v1/gonum/mat/triangular.go index 743fd38c..0e37fb01 100644 --- a/vendor/gonum.org/v1/gonum/mat/triangular.go +++ b/vendor/gonum.org/v1/gonum/mat/triangular.go @@ -779,3 +779,54 @@ func (t *TriDense) DoColNonZero(j int, fn func(i, j int, v float64)) { } } } + +// SolveTo solves a triangular system T * X = B or Tᵀ * X = B where T is an n×n +// triangular matrix represented by the receiver and B is a given n×nrhs matrix. +// If T is non-singular, the result will be stored into dst and nil will be +// returned. If T is singular, the contents of dst will be undefined and a +// Condition error will be returned. +// +// If dst is empty, SolveTo will resize it to n×nrhs. If dst is not empty, +// SolveTo will panic if dst is not n×nrhs. +func (t *TriDense) SolveTo(dst *Dense, trans bool, b Matrix) error { + n, nrhs := b.Dims() + if n != t.mat.N { + panic(ErrShape) + } + + dst.reuseAsNonZeroed(n, nrhs) + bU, bTrans := untranspose(b) + if dst == bU { + if bTrans { + work := getDenseWorkspace(n, nrhs, false) + defer putDenseWorkspace(work) + work.Copy(b) + dst.Copy(work) + } + } else { + if rm, ok := bU.(RawMatrixer); ok { + dst.checkOverlap(rm.RawMatrix()) + } + dst.Copy(b) + } + + transT := blas.NoTrans + if trans { + transT = blas.Trans + } + ok := lapack64.Trtrs(transT, t.mat, dst.mat) + if !ok { + return Condition(math.Inf(1)) + } + + work := getFloat64s(3*n, false) + iwork := getInts(n, false) + cond := lapack64.Trcon(CondNorm, t.mat, work, iwork) + putFloat64s(work) + putInts(iwork) + if cond > ConditionTolerance { + return Condition(cond) + } + + return nil +} diff --git a/vendor/gonum.org/v1/gonum/mat/triband.go b/vendor/gonum.org/v1/gonum/mat/triband.go index 9d6b3401..aa0b51d6 100644 --- a/vendor/gonum.org/v1/gonum/mat/triband.go +++ b/vendor/gonum.org/v1/gonum/mat/triband.go @@ -577,13 +577,23 @@ func (t *TriBandDense) SolveTo(dst *Dense, trans bool, b Matrix) error { if n != t.mat.N { panic(ErrShape) } - if b, ok := b.(RawMatrixer); ok && dst != b { - dst.checkOverlap(b.RawMatrix()) - } + dst.reuseAsNonZeroed(n, nrhs) - if dst != b { + bU, bTrans := untranspose(b) + if dst == bU { + if bTrans { + work := getDenseWorkspace(n, nrhs, false) + defer putDenseWorkspace(work) + work.Copy(b) + dst.Copy(work) + } + } else { + if rm, ok := bU.(RawMatrixer); ok { + dst.checkOverlap(rm.RawMatrix()) + } dst.Copy(b) } + var ok bool if trans { ok = lapack64.Tbtrs(blas.Trans, t.mat, dst.mat) diff --git a/vendor/gonum.org/v1/gonum/mat/tridiag.go b/vendor/gonum.org/v1/gonum/mat/tridiag.go index 80a1bd4a..c001d486 100644 --- a/vendor/gonum.org/v1/gonum/mat/tridiag.go +++ b/vendor/gonum.org/v1/gonum/mat/tridiag.go @@ -234,13 +234,23 @@ func (a *Tridiag) SolveTo(dst *Dense, trans bool, b Matrix) error { if n != a.mat.N { panic(ErrShape) } - if b, ok := b.(RawMatrixer); ok && dst != b { - dst.checkOverlap(b.RawMatrix()) - } + dst.reuseAsNonZeroed(n, nrhs) - if dst != b { + bU, bTrans := untranspose(b) + if dst == bU { + if bTrans { + work := getDenseWorkspace(n, nrhs, false) + defer putDenseWorkspace(work) + work.Copy(b) + dst.Copy(work) + } + } else { + if rm, ok := bU.(RawMatrixer); ok { + dst.checkOverlap(rm.RawMatrix()) + } dst.Copy(b) } + var aCopy Tridiag aCopy.CloneFromTridiag(a) var ok bool diff --git a/vendor/gonum.org/v1/gonum/mat/vector.go b/vendor/gonum.org/v1/gonum/mat/vector.go index 2035e809..5c5d3ff7 100644 --- a/vendor/gonum.org/v1/gonum/mat/vector.go +++ b/vendor/gonum.org/v1/gonum/mat/vector.go @@ -837,3 +837,19 @@ func (v *VecDense) RowViewOf(m RawMatrixer, i int) { v.mat.Data = rm.Data[i*rm.Stride : i*rm.Stride+rm.Cols] v.mat.N = rm.Cols } + +// Permute rearranges the elements of the n-vector v in the receiver as +// specified by the permutation p[0],p[1],...,p[n-1] of the integers 0,...,n-1. +// +// If inverse is false, the given permutation is applied: +// +// v[p[i]] is moved to v[i] for i=0,1,...,n-1. +// +// If inverse is true, the inverse permutation is applied: +// +// v[i] is moved to v[p[i]] for i=0,1,...,n-1. +// +// p must have length n, otherwise Permute will panic. +func (v *VecDense) Permute(p []int, inverse bool) { + v.asDense().PermuteRows(p, inverse) +} diff --git a/vendor/modules.txt b/vendor/modules.txt index ff1e1080..7d7ddd4a 100644 --- a/vendor/modules.txt +++ b/vendor/modules.txt @@ -70,7 +70,7 @@ github.com/golang/protobuf/ptypes github.com/golang/protobuf/ptypes/any github.com/golang/protobuf/ptypes/duration github.com/golang/protobuf/ptypes/timestamp -# github.com/google/go-cmp v0.5.9 +# github.com/google/go-cmp v0.6.0 ## explicit; go 1.13 github.com/google/go-cmp/cmp github.com/google/go-cmp/cmp/internal/diff @@ -176,8 +176,8 @@ golang.org/x/text/secure/bidirule golang.org/x/text/transform golang.org/x/text/unicode/bidi golang.org/x/text/unicode/norm -# gonum.org/v1/gonum v0.13.0 -## explicit; go 1.18 +# gonum.org/v1/gonum v0.15.1 +## explicit; go 1.22 gonum.org/v1/gonum/blas gonum.org/v1/gonum/blas/blas64 gonum.org/v1/gonum/blas/cblas128