Repository: gonum/gonum
Pull request #1106 (2019-10-02 07:48:47)

gonum: add release schedule
Please take a look.

Fixes gonum/gonum#1089.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


kortschak: Is there a reason for the quarterly release?

If we are going to do this, I'd like to drop the oldest supported Go version and only cover 2 releases.
sbinet: > Is there a reason for the quarterly release?

no earth-shattering reason, no.
I think the major driver here is to be aligned w/ Go releases.
I'll amend and update the PR for a 2-releases per year release schedule.

> I'd like to drop the oldest supported Go version and only cover 2 releases.

according to the [Go release policy](https://golang.org/doc/devel/release.html#policy), the version support is N-2 *excluded* (I must admit I thought it was N-2 included.)
so I guess we could also piggy-back on that policy.
(I suspect this would mostly only impact (ie: reduce) our travis-ci build times)

one thing I didn't mention in the PR is what we want to do with `go.mod` dependencies.
I think we should `go get -u foo@latest` for each new Gonum release (at least the `golang.org/x/...` ones) to cut on the accumulated cruft.
thoughts?
kortschak: The other reason for doing N and N-1 is that it allows us to take up recent additions to the standard library and language quicker.

I'm not sure of the need to update deps unless there are bugs identified in them; MVS will already pick the more recent version in client code and we don't actually have that many deps (by design).
sbinet: IMHO, updating deps is the same mechanism than doing N and N-1: it makes sure we advance the global ecosystem together towards N,N-1 versions. (and prevents or limits version skew)

for example, in Go-HEP, I get multiple versions of `golang.org/x/...` repos from a bunch of deps (some through gonum and gonum/plot).
kortschak: >for example, in Go-HEP, I get multiple versions of `golang.org/x/...` repos from a bunch of deps (some through gonum and gonum/plot).

What do you mean by this? The compiled code should have only one version of any given dependency.
sbinet: no, you can have many versions of a given module in a given binary.
I myself may be using mod-v0.2.0 while some of my dependencies may be using mod-v0.1.1.
when one looks at the resulting `go.sum` file, one should see all of those.

_e.g.:_ with current `gonum/gonum`:
```
$> cat gonum.org/v1/gonum/go.sum
[...]
golang.org/x/exp v0.0.0-20180321215751-8460e604b9de
golang.org/x/exp v0.0.0-20180807140117-3d87b88a115f
golang.org/x/exp v0.0.0-20190125153040-c74c464bbbf2
golang.org/x/image v0.0.0-20180708004352-c73c2afc3b81
golang.org/x/tools v0.0.0-20180525024113-a5b4c53f6e8b
golang.org/x/tools v0.0.0-20190206041539-40960b6deb8e
[...]
```

3 different versions of `x/exp`, 2 different versions of `x/tools`.
kortschak: That's not what MVS is doing. There are sums for every possible dep, but only one of them will be included in the actual binary.
sbinet: ~~I beg to differ.~~
ok. re-reading (vgo-mvs)[https://research.swtch.com/vgo-mvs], I may have generalized the rule for `M.N` (in `v.M.N.P`) to `N=0` and `P`.

ok. let's not add this clause about upgrading module deps, and stick with "explicitly updating them when a bug is encountered (or a new feature needed)".
kortschak: Yeah, the definitive test is to query an executable. If you run `goversion -m` on some executable built with modules it shows which modules and versions have gone into its build.

Patch:
diff --git a/.travis.yml b/.travis.yml
index 5635bc34..7adf5b02 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -7,7 +7,6 @@ language: go
 go:
  - 1.13.x
  - 1.12.x
- - 1.11.x
  - master
 
 os:
@@ -37,36 +36,22 @@ matrix:
  fast_finish: true
  exclude:
   - os: osx
-    go: 1.13.x
     env: TAGS="-tags bounds"
   - os: osx
-    go: 1.13.x
     env: TAGS="-tags noasm"
   - os: osx
-    go: 1.13.x
     env: TAGS="-tags appengine"
 
-  - os: osx
-    go: 1.12.x
-  - os: osx
-    go: 1.11.x
   - os: osx
     go: master
 
   - os: windows
-    go: 1.13.x
     env: TAGS="-tags bounds"
   - os: windows
-    go: 1.13.x
     env: TAGS="-tags noasm"
   - os: windows
-    go: 1.13.x
     env: TAGS="-tags appengine"
 
-  - os: windows
-    go: 1.12.x
-  - os: windows
-    go: 1.11.x
   - os: windows
     go: master
  allow_failures:
diff --git a/.travis/script.d/deps.sh b/.travis/script.d/deps.sh
index 4eafe9c7..65dd9a3b 100755
--- a/.travis/script.d/deps.sh
+++ b/.travis/script.d/deps.sh
@@ -3,13 +3,10 @@
 set -ex
 
 # Avoid contaminating the go.mod/go.sum files.
+# TODO(kortschak): Remove when golang/go#30515 is resolved
 WORK=$(mktemp -d)
 pushd $WORK
 
-# Work around go1.11 issue where go get will not work outside a module.
-# TODO(kortschak): Remove when go1.11 is no longer supported.
-go mod init fake || true
-
 # Required for format check.
 go get golang.org/x/tools/cmd/goimports
 # Required for imports check.
@@ -20,9 +17,9 @@ go get gonum.org/v1/tools/cmd/check-copyright
 go get golang.org/x/tools/cmd/cover
 go get github.com/mattn/goveralls
 # Required for dot parser checks.
-${TRAVIS_BUILD_DIR}/.travis/script.d/install-gocc.sh 66c61e91b3657c517a6f89d2837d370e61fb9430
+go get github.com/goccmack/gocc@66c61e9
 
 # Clean up.
-# TODO(kortschak): Remove when go1.11 is no longer supported.
+# TODO(kortschak): Remove when golang/go#30515 is resolved.
 popd
 rm -rf $WORK
diff --git a/.travis/script.d/install-gocc.sh b/.travis/script.d/install-gocc.sh
deleted file mode 100755
index 608ed8a1..00000000
--- a/.travis/script.d/install-gocc.sh
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/bin/bash
-
-set -ex
-
-# TODO(kortschak): Replace this with versioned go get when
-# all our supported versions of Go support it.
-mkdir -p $GOPATH/src/github.com/goccmack
-git clone https://github.com/goccmack/gocc.git $GOPATH/src/github.com/goccmack/gocc
-cd $GOPATH/src/github.com/goccmack/gocc
-git checkout $1
-go install
diff --git a/AUTHORS b/AUTHORS
index b560f3fc..132a341d 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -21,6 +21,7 @@ Clayton Northey <clayton.northey@gmail.com>
 Dan Kortschak <dan.kortschak@adelaide.edu.au> <dan@kortschak.io>
 Daniel Fireman <danielfireman@gmail.com>
 Dario Heinisch <dario.heinisch@gmail.com>
+David Kleiven <davidkleiven446@gmail.com>
 David Samborski <bloggingarrow@gmail.com>
 Davor Kapsa <davor.kapsa@gmail.com>
 DeepMind Technologies
@@ -65,7 +66,9 @@ Martin Diz <github@martindiz.com.ar>
 Matthieu Di Mercurio <matthieu.dimercurio@gmail.com>
 Max Halford <maxhalford25@gmail.com>
 MinJae Kwon <k239507@gmail.com>
+Nathan Edwards <etaoinshrdluwho@gmail.com>
 Nick Potts <nick@the-potts.com>
+Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
 Or Rikon <rikonor@gmail.com>
 Pontus Melke <pontusmelke@gmail.com>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index efb56ec7..e90db65c 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -29,6 +29,7 @@ Clayton Northey <clayton.northey@gmail.com>
 Dan Kortschak <dan.kortschak@adelaide.edu.au> <dan@kortschak.io>
 Daniel Fireman <danielfireman@gmail.com>
 Dario Heinisch <dario.heinisch@gmail.com>
+David Kleiven <davidkleiven446@gmail.com>
 David Samborski <bloggingarrow@gmail.com>
 Davor Kapsa <davor.kapsa@gmail.com>
 Delaney Gillilan <delaneygillilan@gmail.com>
@@ -71,7 +72,9 @@ Martin Diz <github@martindiz.com.ar>
 Matthieu Di Mercurio <matthieu.dimercurio@gmail.com>
 Max Halford <maxhalford25@gmail.com>
 MinJae Kwon <k239507@gmail.com>
+Nathan Edwards <etaoinshrdluwho@gmail.com>
 Nick Potts <nick@the-potts.com>
+Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
 Or Rikon <rikonor@gmail.com>
 Pontus Melke <pontusmelke@gmail.com>
diff --git a/README.md b/README.md
index 81c25b2f..8a5e0b05 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 ## Installation
 
-The core packages of the gonum suite are written in pure Go with some assembly.
+The core packages of the Gonum suite are written in pure Go with some assembly.
 Installation is done using `go get`.
 ```
 go get -u gonum.org/v1/gonum/...
@@ -10,13 +10,23 @@ go get -u gonum.org/v1/gonum/...
 
 ## Supported Go versions
 
-Gonum supports and tests on the three most recent minor versions of Go on [Linux](https://github.com/gonum/gonum/blob/master/.travis.yml#L6-L11) and [Windows](https://github.com/gonum/gonum/blob/master/appveyor.yml#L13-L18).
+Gonum supports and tests on the [two most recent Go releases](https://github.com/gonum/gonum/blob/master/.travis.yml#L6-L10) on Linux, macOS and Windows.
 
+## Release schedule
+
+The Gonum modules are released on a six-month release schedule, aligned with the Go releases.
+_i.e.:_ when `Go-1.x` is released, `Gonum-v0.n.0` is released around the same time.
+Six months after, `Go-1.x+1` is released, and `Gonum-v0.n+1.0` as well.
+
+The release schedule, based on the current Go release schedule is thus:
+
+- `Gonum-v0.n.0`: February
+- `Gonum-v0.n+1.0`: August
 
 ## Build tags
 
-The gonum packages use a variety of build tags to set non-standard build conditions.
-Building gonum applications will work without knowing how to use these tags, but they can be used during testing and to control the use of assembly and CGO code.
+The Gonum packages use a variety of build tags to set non-standard build conditions.
+Building Gonum applications will work without knowing how to use these tags, but they can be used during testing and to control the use of assembly and CGO code.
 
 The current list of non-internal tags is as follows:
 
@@ -38,4 +48,4 @@ https://groups.google.com/forum/#!forum/gonum-dev
 
 Original code is licensed under the Gonum License found in the LICENSE file. Portions of the code are subject to the additional licenses found in THIRD_PARTY_LICENSES. All third party code is licensed either under a BSD or MIT license.
 
-Code in graph/formats/dot is dual licensed [Public Domain Dedication](https://creativecommons.org/publicdomain/zero/1.0/) and Gonum License, and users are free to choose the license which suits their needs for this code.
\ No newline at end of file
+Code in graph/formats/dot is dual licensed [Public Domain Dedication](https://creativecommons.org/publicdomain/zero/1.0/) and Gonum License, and users are free to choose the license which suits their needs for this code.
diff --git a/appveyor.yml b/appveyor.yml
index 94009f7d..003bce1b 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -13,7 +13,6 @@ environment:
   # Do not move these lines; they are referred to by README.md.
   # Versions of go that are explicitly supported by gonum.
   matrix:
-    - GOROOT: 'c:\go111'
     - GOROOT: 'c:\go112'
     #- GOROOT: 'c:\go113'
   GOPATH: c:\gopath
diff --git a/graph/encoding/dot/encode.go b/graph/encoding/dot/encode.go
index 10fc7d95..3e3955c5 100644
--- a/graph/encoding/dot/encode.go
+++ b/graph/encoding/dot/encode.go
@@ -135,8 +135,6 @@ type printer struct {
 	prefix string
 	indent string
 	depth  int
-
-	err error
 }
 
 type edge struct {
@@ -291,7 +289,7 @@ func (p *simpleGraphPrinter) print(g graph.Graph, name string, needsIndent, isSu
 	return nil
 }
 
-func (p *printer) printFrontMatter(name string, needsIndent, isSubgraph, isDirected, isStrict bool) error {
+func (p *printer) printFrontMatter(name string, needsIndent, isSubgraph, isDirected, isStrict bool) {
 	p.buf.WriteString(p.prefix)
 	if needsIndent {
 		for i := 0; i < p.depth; i++ {
@@ -316,7 +314,6 @@ func (p *printer) printFrontMatter(name string, needsIndent, isSubgraph, isDirec
 	}
 
 	p.openBlock(" {")
-	return nil
 }
 
 func (p *printer) writeNode(n graph.Node) {
diff --git a/graph/ex/fdpclust/gn.go b/graph/ex/fdpclust/gn.go
deleted file mode 100644
index 4df32b0d..00000000
--- a/graph/ex/fdpclust/gn.go
+++ /dev/null
@@ -1,282 +0,0 @@
-// Copyright ©2014 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 main
-
-import (
-	"gonum.org/v1/gonum/graph"
-	"gonum.org/v1/gonum/graph/iterator"
-	"gonum.org/v1/gonum/graph/simple"
-)
-
-type GraphNode struct {
-	id        int64
-	neighbors []graph.Node
-	roots     []*GraphNode
-}
-
-func (g *GraphNode) Node(id int64) graph.Node {
-	if id == g.id {
-		return g
-	}
-
-	visited := map[int64]struct{}{g.id: {}}
-	for _, root := range g.roots {
-		if root.ID() == id {
-			return root
-		}
-
-		if root.has(id, visited) {
-			return root
-		}
-	}
-
-	for _, neigh := range g.neighbors {
-		if neigh.ID() == id {
-			return neigh
-		}
-
-		if gn, ok := neigh.(*GraphNode); ok {
-			if gn.has(id, visited) {
-				return gn
-			}
-		}
-	}
-
-	return nil
-}
-
-func (g *GraphNode) has(id int64, visited map[int64]struct{}) bool {
-	for _, root := range g.roots {
-		if _, ok := visited[root.ID()]; ok {
-			continue
-		}
-
-		visited[root.ID()] = struct{}{}
-		if root.ID() == id {
-			return true
-		}
-
-		if root.has(id, visited) {
-			return true
-		}
-
-	}
-
-	for _, neigh := range g.neighbors {
-		if _, ok := visited[neigh.ID()]; ok {
-			continue
-		}
-
-		visited[neigh.ID()] = struct{}{}
-		if neigh.ID() == id {
-			return true
-		}
-
-		if gn, ok := neigh.(*GraphNode); ok {
-			if gn.has(id, visited) {
-				return true
-			}
-		}
-	}
-
-	return false
-}
-
-func (g *GraphNode) Nodes() graph.Nodes {
-	toReturn := []graph.Node{g}
-	visited := map[int64]struct{}{g.id: {}}
-
-	for _, root := range g.roots {
-		toReturn = append(toReturn, root)
-		visited[root.ID()] = struct{}{}
-
-		toReturn = root.nodes(toReturn, visited)
-	}
-
-	for _, neigh := range g.neighbors {
-		toReturn = append(toReturn, neigh)
-		visited[neigh.ID()] = struct{}{}
-
-		if gn, ok := neigh.(*GraphNode); ok {
-			toReturn = gn.nodes(toReturn, visited)
-		}
-	}
-
-	return iterator.NewOrderedNodes(toReturn)
-}
-
-func (g *GraphNode) nodes(list []graph.Node, visited map[int64]struct{}) []graph.Node {
-	for _, root := range g.roots {
-		if _, ok := visited[root.ID()]; ok {
-			continue
-		}
-		visited[root.ID()] = struct{}{}
-		list = append(list, graph.Node(root))
-
-		list = root.nodes(list, visited)
-	}
-
-	for _, neigh := range g.neighbors {
-		if _, ok := visited[neigh.ID()]; ok {
-			continue
-		}
-
-		list = append(list, neigh)
-		if gn, ok := neigh.(*GraphNode); ok {
-			list = gn.nodes(list, visited)
-		}
-	}
-
-	return list
-}
-
-func (g *GraphNode) From(id int64) graph.Nodes {
-	if id == g.ID() {
-		return iterator.NewOrderedNodes(g.neighbors)
-	}
-
-	visited := map[int64]struct{}{g.id: {}}
-	for _, root := range g.roots {
-		visited[root.ID()] = struct{}{}
-
-		if result := root.findNeighbors(id, visited); result != nil {
-			return iterator.NewOrderedNodes(result)
-		}
-	}
-
-	for _, neigh := range g.neighbors {
-		visited[neigh.ID()] = struct{}{}
-
-		if gn, ok := neigh.(*GraphNode); ok {
-			if result := gn.findNeighbors(id, visited); result != nil {
-				return iterator.NewOrderedNodes(result)
-			}
-		}
-	}
-
-	return nil
-}
-
-func (g *GraphNode) findNeighbors(id int64, visited map[int64]struct{}) []graph.Node {
-	if id == g.ID() {
-		return g.neighbors
-	}
-
-	for _, root := range g.roots {
-		if _, ok := visited[root.ID()]; ok {
-			continue
-		}
-		visited[root.ID()] = struct{}{}
-
-		if result := root.findNeighbors(id, visited); result != nil {
-			return result
-		}
-	}
-
-	for _, neigh := range g.neighbors {
-		if _, ok := visited[neigh.ID()]; ok {
-			continue
-		}
-		visited[neigh.ID()] = struct{}{}
-
-		if gn, ok := neigh.(*GraphNode); ok {
-			if result := gn.findNeighbors(id, visited); result != nil {
-				return result
-			}
-		}
-	}
-
-	return nil
-}
-
-func (g *GraphNode) HasEdgeBetween(uid, vid int64) bool {
-	return g.EdgeBetween(uid, vid) != nil
-}
-
-func (g *GraphNode) Edge(uid, vid int64) graph.Edge {
-	return g.EdgeBetween(uid, vid)
-}
-
-func (g *GraphNode) EdgeBetween(uid, vid int64) graph.Edge {
-	if uid == g.id || vid == g.id {
-		for _, neigh := range g.neighbors {
-			if neigh.ID() == uid || neigh.ID() == vid {
-				return simple.Edge{F: g, T: neigh}
-			}
-		}
-		return nil
-	}
-
-	visited := map[int64]struct{}{g.id: {}}
-	for _, root := range g.roots {
-		visited[root.ID()] = struct{}{}
-		if result := root.edgeBetween(uid, vid, visited); result != nil {
-			return result
-		}
-	}
-
-	for _, neigh := range g.neighbors {
-		visited[neigh.ID()] = struct{}{}
-		if gn, ok := neigh.(*GraphNode); ok {
-			if result := gn.edgeBetween(uid, vid, visited); result != nil {
-				return result
-			}
-		}
-	}
-
-	return nil
-}
-
-func (g *GraphNode) edgeBetween(uid, vid int64, visited map[int64]struct{}) graph.Edge {
-	if uid == g.id || vid == g.id {
-		for _, neigh := range g.neighbors {
-			if neigh.ID() == uid || neigh.ID() == vid {
-				return simple.Edge{F: g, T: neigh}
-			}
-		}
-		return nil
-	}
-
-	for _, root := range g.roots {
-		if _, ok := visited[root.ID()]; ok {
-			continue
-		}
-		visited[root.ID()] = struct{}{}
-		if result := root.edgeBetween(uid, vid, visited); result != nil {
-			return result
-		}
-	}
-
-	for _, neigh := range g.neighbors {
-		if _, ok := visited[neigh.ID()]; ok {
-			continue
-		}
-
-		visited[neigh.ID()] = struct{}{}
-		if gn, ok := neigh.(*GraphNode); ok {
-			if result := gn.edgeBetween(uid, vid, visited); result != nil {
-				return result
-			}
-		}
-	}
-
-	return nil
-}
-
-func (g *GraphNode) ID() int64 {
-	return g.id
-}
-
-func (g *GraphNode) AddNeighbor(n *GraphNode) {
-	g.neighbors = append(g.neighbors, graph.Node(n))
-}
-
-func (g *GraphNode) AddRoot(n *GraphNode) {
-	g.roots = append(g.roots, n)
-}
-
-func NewGraphNode(id int64) *GraphNode {
-	return &GraphNode{id: id}
-}
diff --git a/graph/ex/fdpclust/main.go b/graph/ex/fdpclust/main.go
deleted file mode 100644
index a1613e80..00000000
--- a/graph/ex/fdpclust/main.go
+++ /dev/null
@@ -1,79 +0,0 @@
-// Copyright ©2014 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 main
-
-import (
-	"fmt"
-
-	"gonum.org/v1/gonum/graph"
-	"gonum.org/v1/gonum/graph/topo"
-)
-
-func main() {
-	// graph G {
-	G := NewGraphNode(0)
-	// e
-	e := NewGraphNode(1)
-
-	// subgraph clusterA {
-	clusterA := NewGraphNode(2)
-
-	// a -- b
-	a := NewGraphNode(3)
-	b := NewGraphNode(4)
-	a.AddNeighbor(b)
-	b.AddNeighbor(a)
-	clusterA.AddRoot(a)
-	clusterA.AddRoot(b)
-
-	// subgraph clusterC {
-	clusterC := NewGraphNode(5)
-	// C -- D
-	C := NewGraphNode(6)
-	D := NewGraphNode(7)
-	C.AddNeighbor(D)
-	D.AddNeighbor(C)
-
-	clusterC.AddRoot(C)
-	clusterC.AddRoot(D)
-	// }
-	clusterA.AddRoot(clusterC)
-	// }
-
-	// subgraph clusterB {
-	clusterB := NewGraphNode(8)
-
-	// d -- f
-	d := NewGraphNode(9)
-	f := NewGraphNode(10)
-	d.AddNeighbor(f)
-	f.AddNeighbor(d)
-	clusterB.AddRoot(d)
-	clusterB.AddRoot(f)
-	// }
-
-	// d -- D
-	d.AddNeighbor(D)
-	D.AddNeighbor(d)
-
-	// e -- clusterB
-	e.AddNeighbor(clusterB)
-	clusterB.AddNeighbor(e)
-
-	// clusterC -- clusterB
-	clusterC.AddNeighbor(clusterB)
-	clusterB.AddNeighbor(clusterC)
-
-	G.AddRoot(e)
-	G.AddRoot(clusterA)
-	G.AddRoot(clusterB)
-	// }
-
-	if !topo.IsPathIn(G, []graph.Node{C, D, d, f}) {
-		fmt.Println("Not working!")
-	} else {
-		fmt.Println("Working!")
-	}
-}
diff --git a/graph/iterator/nodes.go b/graph/iterator/nodes.go
index 952dd770..f32ec121 100644
--- a/graph/iterator/nodes.go
+++ b/graph/iterator/nodes.go
@@ -4,7 +4,11 @@
 
 package iterator
 
-import "gonum.org/v1/gonum/graph"
+import (
+	"reflect"
+
+	"gonum.org/v1/gonum/graph"
+)
 
 // OrderedNodes implements the graph.Nodes and graph.NodeSlicer interfaces.
 // The iteration order of OrderedNodes is the order of nodes passed to
@@ -123,3 +127,54 @@ func (n *ImplicitNodes) NodeSlice() []graph.Node {
 	}
 	return nodes
 }
+
+// Nodes implements the graph.Nodes interfaces.
+// The iteration order of Nodes is randomized.
+type Nodes struct {
+	nodes reflect.Value
+	iter  *reflect.MapIter
+	pos   int
+	curr  graph.Node
+}
+
+// NewNodes returns a Nodes initialized with the provided nodes, a
+// map of node IDs to graph.Nodes. No check is made that the keys
+// match the graph.Node IDs, and the map keys are not used.
+//
+// Behavior of the Nodes is unspecified if nodes is mutated after
+// the call the NewNodes.
+func NewNodes(nodes map[int64]graph.Node) *Nodes {
+	rv := reflect.ValueOf(nodes)
+	return &Nodes{nodes: rv, iter: rv.MapRange()}
+}
+
+// Len returns the remaining number of nodes to be iterated over.
+func (n *Nodes) Len() int {
+	return n.nodes.Len() - n.pos
+}
+
+// Next returns whether the next call of Node will return a valid node.
+func (n *Nodes) Next() bool {
+	if n.pos >= n.nodes.Len() {
+		return false
+	}
+	ok := n.iter.Next()
+	if ok {
+		n.pos++
+		n.curr = n.iter.Value().Interface().(graph.Node)
+	}
+	return ok
+}
+
+// Node returns the current node of the iterator. Next must have been
+// called prior to a call to Node.
+func (n *Nodes) Node() graph.Node {
+	return n.curr
+}
+
+// Reset returns the iterator to its initial state.
+func (n *Nodes) Reset() {
+	n.curr = nil
+	n.pos = 0
+	n.iter = n.nodes.MapRange()
+}
diff --git a/graph/iterator/nodes_unordered.go b/graph/iterator/nodes_unordered.go
deleted file mode 100644
index a6232f9c..00000000
--- a/graph/iterator/nodes_unordered.go
+++ /dev/null
@@ -1,66 +0,0 @@
-// Copyright ©2018 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.
-
-// +build go1.12
-
-// TODO(kortschak): Merge this into nodes.go when go1.11 is no longer supported.
-
-package iterator
-
-import (
-	"reflect"
-
-	"gonum.org/v1/gonum/graph"
-)
-
-// Nodes implements the graph.Nodes interfaces.
-// The iteration order of Nodes is randomized.
-type Nodes struct {
-	nodes reflect.Value
-	iter  *reflect.MapIter
-	pos   int
-	curr  graph.Node
-}
-
-// NewNodes returns a Nodes initialized with the provided nodes, a
-// map of node IDs to graph.Nodes. No check is made that the keys
-// match the graph.Node IDs, and the map keys are not used.
-//
-// Behavior of the Nodes is unspecified if nodes is mutated after
-// the call the NewNodes.
-func NewNodes(nodes map[int64]graph.Node) *Nodes {
-	rv := reflect.ValueOf(nodes)
-	return &Nodes{nodes: rv, iter: rv.MapRange()}
-}
-
-// Len returns the remaining number of nodes to be iterated over.
-func (n *Nodes) Len() int {
-	return n.nodes.Len() - n.pos
-}
-
-// Next returns whether the next call of Node will return a valid node.
-func (n *Nodes) Next() bool {
-	if n.pos >= n.nodes.Len() {
-		return false
-	}
-	ok := n.iter.Next()
-	if ok {
-		n.pos++
-		n.curr = n.iter.Value().Interface().(graph.Node)
-	}
-	return ok
-}
-
-// Node returns the current node of the iterator. Next must have been
-// called prior to a call to Node.
-func (n *Nodes) Node() graph.Node {
-	return n.curr
-}
-
-// Reset returns the iterator to its initial state.
-func (n *Nodes) Reset() {
-	n.curr = nil
-	n.pos = 0
-	n.iter = n.nodes.MapRange()
-}
diff --git a/graph/layout/isomap.go b/graph/layout/isomap.go
index 4ce9cf84..879f6c37 100644
--- a/graph/layout/isomap.go
+++ b/graph/layout/isomap.go
@@ -53,11 +53,11 @@ func isomap(g graph.Graph, nodes []graph.Node, dims int) *mat.Dense {
 			dist.SetSym(i, j, p.Weight(u.ID(), v.ID()))
 		}
 	}
-
-	k, v, _ := mds.TorgersonScaling(nil, nil, dist)
+	var v mat.Dense
+	k, _ := mds.TorgersonScaling(&v, nil, dist)
 	if k < dims {
 		return nil
 	}
 
-	return v
+	return &v
 }
diff --git a/graph/path/a_star.go b/graph/path/a_star.go
index e12d655b..c502abd3 100644
--- a/graph/path/a_star.go
+++ b/graph/path/a_star.go
@@ -9,6 +9,7 @@ import (
 
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/set"
+	"gonum.org/v1/gonum/graph/traverse"
 )
 
 // AStar finds the A*-shortest path from s to t in g using the heuristic h. The path and
@@ -23,9 +24,11 @@ import (
 // If h is nil, AStar will use the g.HeuristicCost method if g implements HeuristicCoster,
 // falling back to NullHeuristic otherwise. If the graph does not implement Weighted,
 // UniformCost is used. AStar will panic if g has an A*-reachable negative edge weight.
-func AStar(s, t graph.Node, g graph.Graph, h Heuristic) (path Shortest, expanded int) {
-	if g.Node(s.ID()) == nil || g.Node(t.ID()) == nil {
-		return Shortest{from: s}, 0
+func AStar(s, t graph.Node, g traverse.Graph, h Heuristic) (path Shortest, expanded int) {
+	if g, ok := g.(graph.Graph); ok {
+		if g.Node(s.ID()) == nil || g.Node(t.ID()) == nil {
+			return Shortest{from: s}, 0
+		}
 	}
 	var weight Weighting
 	if wg, ok := g.(Weighted); ok {
@@ -41,7 +44,7 @@ func AStar(s, t graph.Node, g graph.Graph, h Heuristic) (path Shortest, expanded
 		}
 	}
 
-	path = newShortestFrom(s, graph.NodesOf(g.Nodes()))
+	path = newShortestFrom(s, []graph.Node{s, t})
 	tid := t.ID()
 
 	visited := make(set.Int64s)
@@ -64,14 +67,17 @@ func AStar(s, t graph.Node, g graph.Graph, h Heuristic) (path Shortest, expanded
 			if visited.Has(vid) {
 				continue
 			}
-			j := path.indexOf[vid]
+			j, ok := path.indexOf[vid]
+			if !ok {
+				j = path.add(v)
+			}
 
 			w, ok := weight(u.node.ID(), vid)
 			if !ok {
-				panic("A*: unexpected invalid weight")
+				panic("path: A* unexpected invalid weight")
 			}
 			if w < 0 {
-				panic("A*: negative edge weight")
+				panic("path: A* negative edge weight")
 			}
 			g := u.gscore + w
 			if n, ok := open.node(vid); !ok {
diff --git a/graph/path/dynamic/dstarlite.go b/graph/path/dynamic/dstarlite.go
index 3dc0acce..f388b698 100644
--- a/graph/path/dynamic/dstarlite.go
+++ b/graph/path/dynamic/dstarlite.go
@@ -406,14 +406,18 @@ queue:
 // key is a D* Lite priority queue key.
 type key [2]float64
 
+// badKey is a poisoned key. Testing for a bad key uses NaN inequality.
 var badKey = key{math.NaN(), math.NaN()}
 
+//nolint:staticcheck
+func (k key) isBadKey() bool { return k != k }
+
 // less returns whether k is less than other. From ISBN:0-262-51129-0 pp476-483:
 //
 //  k ≤ k' iff k₁ < k'₁ OR (k₁ == k'₁ AND k₂ ≤ k'₂)
 //
 func (k key) less(other key) bool {
-	if k != k || other != other {
+	if k.isBadKey() || other.isBadKey() {
 		panic("D* Lite: poisoned key")
 	}
 	return k[0] < other[0] || (k[0] == other[0] && k[1] < other[1])
diff --git a/integrate/romberg.go b/integrate/romberg.go
new file mode 100644
index 00000000..25e45d0a
--- /dev/null
+++ b/integrate/romberg.go
@@ -0,0 +1,65 @@
+// Copyright ©2019 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 integrate
+
+import (
+	"math"
+	"math/bits"
+)
+
+// Romberg returns an approximate value of the integral
+//  \int_a^b f(x)dx
+// computed using the Romberg's method. The function f is given
+// as a slice of equally-spaced samples, that is,
+//  f[i] = f(a + i*dx)
+// and dx is the spacing between the samples.
+//
+// The length of f must be 2^k + 1, where k is a positive integer,
+// and dx must be positive.
+//
+// See https://en.wikipedia.org/wiki/Romberg%27s_method for a description of
+// the algorithm.
+func Romberg(f []float64, dx float64) float64 {
+	if len(f) < 3 {
+		panic("integral: invalid slice length: must be at least 3")
+	}
+
+	if dx <= 0 {
+		panic("integral: invalid spacing: must be larger than 0")
+	}
+
+	n := len(f) - 1
+	k := bits.Len(uint(n - 1))
+
+	if len(f) != 1<<uint(k)+1 {
+		panic("integral: invalid slice length: must be 2^k + 1")
+	}
+
+	work := make([]float64, 2*(k+1))
+	prev := work[:k+1]
+	curr := work[k+1:]
+
+	h := dx * float64(n)
+	prev[0] = (f[0] + f[n]) * 0.5 * h
+
+	step := n
+	for i := 1; i <= k; i++ {
+		h /= 2
+		step /= 2
+		var estimate float64
+		for j := 0; j < n/2; j += step {
+			estimate += f[2*j+step]
+		}
+
+		curr[0] = estimate*h + 0.5*prev[0]
+		for j := 1; j <= i; j++ {
+			factor := math.Pow(4, float64(j))
+			curr[j] = (factor*curr[j-1] - prev[j-1]) / (factor - 1)
+		}
+
+		prev, curr = curr, prev
+	}
+	return prev[k]
+}
diff --git a/integrate/simpsons.go b/integrate/simpsons.go
new file mode 100644
index 00000000..35a9e2ae
--- /dev/null
+++ b/integrate/simpsons.go
@@ -0,0 +1,93 @@
+// Copyright ©2019 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 integrate
+
+import "sort"
+
+// Simpsons returns an approximate value of the integral
+//  \int_a^b f(x)dx
+// computed using the Simpsons's method. The function f is given as a slice of
+// samples evaluated at locations in x, that is,
+//  f[i] = f(x[i]), x[0] = a, x[len(x)-1] = b
+// The slice x must be sorted in strictly increasing order. x and f must be of
+// equal length and the length must be at least 3.
+//
+// See https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
+// for more information.
+func Simpsons(x, f []float64) float64 {
+	n := len(x)
+	switch {
+	case len(f) != n:
+		panic("integrate: slice length mismatch")
+	case n < 3:
+		panic("integrate: input data too small")
+	case !sort.Float64sAreSorted(x):
+		panic("integrate: must be sorted")
+	}
+
+	// Simpson's method approximates the integral of a function f by subdividing
+	// the interval of integration into segments and applying a method that fits
+	// a polynomial to each subinterval.
+	//
+	// This implementation makes use of the following composite Simpson's method
+	// to estimate \int_a^b f(x) dx where a, b are x at x[0] and x[len(x)-1]
+	// respectively:
+	//  \sum_{i=1}^{N} {a_0}_i * f_{i-1} + {a_1}_i * f_{i} + {a_2}_i * f_{i+1}
+	// where N is the count of subintervals and {a_0}_i, {a_1}_i, and {a_2}_i
+	// are constants at index i given by:
+	//  {a_0}_i = 2 * h^3_i - h^3_{i+1} + 3 * h_{i+1} * h^2_i / 6 * h_i * (h_i + h_{i+1})
+	//  {a_1}_i = h^3_i + h^3_{i+1} + 3 * h_i * h_{i+1} * (h_i + h_{i+1}) / 6 * h_i * h_{i+1}
+	//  {a_2}_i = -h^3_i + 2 * h^3_{i+1} + 3 * h_i * h^2_{i+1} / 6 * h_{i+1} * (h_i + h_{i+1})
+	// where h_k is the difference x[k] - x[k-1].
+	//
+	// The formula above approximates the integral of f when N is even. If the
+	// number of subintervals is odd, the subintervals i=0,...,n-2 are given by
+	// the above formula and the approximation over the remaining subintervals
+	// is given by:
+	//  {a_0}_i * f_{N-2} + {a_1}_i * f_{N-1} + {a_2}_i * f_N
+	// where the coefficients are:
+	//  {a_0}_i = -1 * h^3_{N-1} / 6 * h_{N-2} * (h_{N-2} + h_{N-1})
+	//  {a_1}_i = h^2_{N-1} + 3 * h^3_{N-1} * h_{N-2} / 6 * h_{N-2}
+	//  {a_2}_i = 2 * h^2_{N-1} + 3 * h^2_{N-1} * h_{N-2} / 6 * (h_{N-2} + h_{N-1})
+
+	var integral float64
+	for i := 1; i < n-1; i += 2 {
+		h0 := x[i] - x[i-1]
+		h1 := x[i+1] - x[i]
+		if h0 == 0 || h1 == 0 {
+			panic("integrate: repeated abscissa")
+		}
+		h0p2 := h0 * h0
+		h0p3 := h0 * h0 * h0
+		h1p2 := h1 * h1
+		h1p3 := h1 * h1 * h1
+		hph := h0 + h1
+		a0 := (2*h0p3 - h1p3 + 3*h1*h0p2) / (6 * h0 * hph)
+		a1 := (h0p3 + h1p3 + 3*h0*h1*hph) / (6 * h0 * h1)
+		a2 := (-h0p3 + 2*h1p3 + 3*h0*h1p2) / (6 * h1 * hph)
+		integral += a0 * f[i-1]
+		integral += a1 * f[i]
+		integral += a2 * f[i+1]
+	}
+
+	if n%2 == 0 {
+		h0 := x[n-2] - x[n-3]
+		h1 := x[n-1] - x[n-2]
+		if h0 == 0 || h1 == 0 {
+			panic("integrate: repeated abscissa")
+		}
+		h1p2 := h1 * h1
+		h1p3 := h1 * h1 * h1
+		hph := h0 + h1
+		a0 := -1 * h1p3 / (6 * h0 * hph)
+		a1 := (h1p2 + 3*h0*h1) / (6 * h0)
+		a2 := (2*h1p2 + 3*h0*h1) / (6 * hph)
+		integral += a0 * f[n-3]
+		integral += a1 * f[n-2]
+		integral += a2 * f[n-1]
+	}
+
+	return integral
+}
diff --git a/integrate/trapezoidal.go b/integrate/trapezoidal.go
index a16096d0..1b8e4a13 100644
--- a/integrate/trapezoidal.go
+++ b/integrate/trapezoidal.go
@@ -6,33 +6,34 @@ package integrate
 
 import "sort"
 
-// Trapezoidal estimates the integral of a function f
+// Trapezoidal returns an approximate value of the integral
 //  \int_a^b f(x) dx
-// from a set of evaluations of the function using the trapezoidal rule.
-// The trapezoidal rule makes piecewise linear approximations to the function,
-// and estimates
+// computed using the trapezoidal rule. The function f is given as a slice of
+// samples evaluated at locations in x, that is,
+//  f[i] = f(x[i]), x[0] = a, x[len(x)-1] = b
+// The slice x must be sorted in strictly increasing order. x and f must be of
+// equal length and the length must be at least 2.
+//
+// The trapezoidal rule approximates f by a piecewise linear function and
+// estimates
 //  \int_x[i]^x[i+1] f(x) dx
 // as
 //  (x[i+1] - x[i]) * (f[i] + f[i+1])/2
-// where f[i] is the value of the function at x[i].
 // More details on the trapezoidal rule can be found at:
 // https://en.wikipedia.org/wiki/Trapezoidal_rule
-//
-// The (x,f) input data points must be sorted along x. stat.SortWeighted can
-// be used to do this. The x and f slices must be of equal length and have
-// length greater than one.
 func Trapezoidal(x, f []float64) float64 {
+	n := len(x)
 	switch {
-	case len(x) != len(f):
+	case len(f) != n:
 		panic("integrate: slice length mismatch")
-	case len(x) < 2:
+	case n < 2:
 		panic("integrate: input data too small")
 	case !sort.Float64sAreSorted(x):
 		panic("integrate: input must be sorted")
 	}
 
 	integral := 0.0
-	for i := 0; i < len(x)-1; i++ {
+	for i := 0; i < n-1; i++ {
 		integral += 0.5 * (x[i+1] - x[i]) * (f[i+1] + f[i])
 	}
 
diff --git a/mathext/airy.go b/mathext/airy.go
index c94c5d84..7b8865f1 100644
--- a/mathext/airy.go
+++ b/mathext/airy.go
@@ -17,7 +17,7 @@ func AiryAi(z complex128) complex128 {
 	// documentation for the exact behavior.
 	id := 0
 	kode := 1
-	air, aii, _ := amos.Zairy(real(z), imag(z), id, kode)
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
 	return complex(air, aii)
 }
 
@@ -32,6 +32,6 @@ func AiryAiDeriv(z complex128) complex128 {
 	// documentation for the exact behavior.
 	id := 1
 	kode := 1
-	air, aii, _ := amos.Zairy(real(z), imag(z), id, kode)
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
 	return complex(air, aii)
 }
diff --git a/mathext/internal/amos/amos.go b/mathext/internal/amos/amos.go
index 35648ceb..57f2bdec 100644
--- a/mathext/internal/amos/amos.go
+++ b/mathext/internal/amos/amos.go
@@ -42,7 +42,7 @@ func max(a, b int) int {
 	return b
 }
 
-func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
+func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ, IERR int) {
 	// zairy is adapted from the original Netlib code by Donald Amos.
 	// http://www.netlib.no/netlib/amos/zairy.f
 
@@ -179,7 +179,7 @@ func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
 		DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
 		S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
 		ZEROR, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB float64
-	var IERR, IFLAG, K, K1, K2, MR, NN int
+	var IFLAG, K, K1, K2, MR, NN int
 	var tmp complex128
 
 	// Extra element for padding.
@@ -420,7 +420,7 @@ OneHundred:
 	if ZI < 0.0e0 {
 		MR = -1
 	}
-	ZTAR, ZTAI, FNU, KODE, MR, _, CYR, CYI, NN, RL, TOL, ELIM, ALIM = Zacai(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, ELIM, ALIM)
+	_, _, _, _, _, _, CYR, CYI, NN, _, _, _, _ = Zacai(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, RL, TOL, ELIM, ALIM)
 	if NN < 0 {
 		goto TwoEighty
 	}
@@ -443,7 +443,7 @@ OneTen:
 		goto TwoTen
 	}
 OneTwenty:
-	ZTAR, ZTAI, FNU, KODE, _, CYR, CYI, NZ, TOL, ELIM, ALIM = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM)
+	_, _, _, _, _, CYR, CYI, NZ, _, _, _ = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, TOL, ELIM, ALIM)
 
 OneThirty:
 	S1R = CYR[1] * COEF
@@ -532,7 +532,7 @@ TwoSixty:
 }
 
 // sbknu computes the k bessel function in the right half z plane.
-func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM, ALIM float64) (ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout, ELIMout, ALIMout float64) {
+func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, TOL, ELIM, ALIM float64) (ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, TOLout, ELIMout, ALIMout float64) {
 	/* Old dimension comment.
 		DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
 	     * CYI(2)
@@ -591,7 +591,6 @@ func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM
 	BRY[1] = 1.0e+3 * dmach[1] / TOL
 	BRY[2] = 1.0e0 / BRY[1]
 	BRY[3] = dmach[2]
-	NZ = 0
 	IFLAG = 0
 	KODED = KODE
 	RCAZ = 1.0e0 / CAZ
@@ -1183,7 +1182,7 @@ TwoSeventy:
 	YI[2] = S2I
 TwoEighty:
 	ASCLE = BRY[1]
-	ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM = Zkscl(ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM)
+	_, _, FNU, N, YR, YI, NZ, RZR, RZI, _, TOL, ELIM = Zkscl(ZDR, ZDI, FNU, N, YR, YI, RZR, RZI, ASCLE, TOL, ELIM)
 	INU = N - NZ
 	if INU <= 0 {
 		return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM
@@ -1213,7 +1212,6 @@ TwoNinety:
 
 	// SCALE BY math.Exp(Z), IFLAG = 1 CASES.
 
-	KODED = 2
 	IFLAG = 1
 	KFLAG = 2
 	goto OneTwenty
@@ -1234,8 +1232,8 @@ ThreeTen:
 // SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
 // ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
 // return WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
-func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, NZ int, RZR, RZI, ASCLE, TOL, ELIM float64) (
-	ZRRout, ZRIout, FNUout float64, Nout int, YRout, YIout []float64, NZout int, RZRout, RZIout, ASCLEout, TOLout, ELIMout float64) {
+func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, RZR, RZI, ASCLE, TOL, ELIM float64) (
+	ZRRout, ZRIout, FNUout float64, Nout int, YRout, YIout []float64, NZ int, RZRout, RZIout, ASCLEout, TOLout, ELIMout float64) {
 	var ACS, AS, CKI, CKR, CSI, CSR, FN, STR, S1I, S1R, S2I,
 		S2R, ZEROI, ZEROR, ZDR, ZDI, CELMR, ELM, HELIM, ALAS float64
 
@@ -1247,7 +1245,6 @@ func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, NZ int, RZR, RZI, ASC
 	// DIMENSION YR(N), YI(N), CYR(2), CYI(2)
 	ZEROR = 0
 	ZEROI = 0
-	NZ = 0
 	IC = 0
 	NN = min(2, N)
 	for I = 1; I <= NN; I++ {
@@ -1413,8 +1410,8 @@ func Zuchk(y complex128, scale, tol float64) int {
 // ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
 // RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT if ZACON
 // IS CALLED FROM ZAIRY.
-func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) (
-	ZRout, ZIout, FNUout float64, KODEout, MRout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) {
+func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, RL, TOL, ELIM, ALIM float64) (
+	ZRout, ZIout, FNUout float64, KODEout, MRout, Nout int, YRout, YIout []float64, NZ int, RLout, TOLout, ELIMout, ALIMout float64) {
 	var ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
 		CSPNI, C1R, C1I, C2R, C2I, DFNU, FMR, PI,
 		SGN, YY, ZNR, ZNI float64
@@ -1427,7 +1424,6 @@ func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, NZ int, RL, T
 	CYI := []float64{math.NaN(), 0, 0}
 
 	PI = math.Pi
-	NZ = 0
 	ZNR = -ZR
 	ZNI = -ZI
 	AZ = cmplx.Abs(complex(ZR, ZI))
@@ -1446,7 +1442,7 @@ Ten:
 	for i, v := range YR {
 		y[i] = complex(v, YI[i])
 	}
-	NW = Zseri(z, FNU, KODE, NN, y[1:], TOL, ELIM, ALIM)
+	Zseri(z, FNU, KODE, NN, y[1:], TOL, ELIM, ALIM)
 	for i, v := range y {
 		YR[i] = real(v)
 		YI[i] = imag(v)
@@ -1457,20 +1453,20 @@ Twenty:
 		goto Thirty
 	}
 	// ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION.
-	ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM = Zasyi(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM)
+	ZNR, ZNI, FNU, KODE, _, YR, YI, NW, RL, TOL, ELIM, ALIM = Zasyi(ZNR, ZNI, FNU, KODE, NN, YR, YI, RL, TOL, ELIM, ALIM)
 	if NW < 0 {
 		goto Eighty
 	}
 	goto Forty
 Thirty:
 	// MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
-	ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL = Zmlri(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
+	ZNR, ZNI, FNU, KODE, _, YR, YI, NW, TOL = Zmlri(ZNR, ZNI, FNU, KODE, NN, YR, YI, TOL)
 	if NW < 0 {
 		goto Eighty
 	}
 Forty:
 	// ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION.
-	ZNR, ZNI, FNU, KODE, _, CYR, CYI, NW, TOL, ELIM, ALIM = Zbknu(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM)
+	ZNR, ZNI, FNU, KODE, _, CYR, CYI, NW, TOL, ELIM, ALIM = Zbknu(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, TOL, ELIM, ALIM)
 	if NW != 0 {
 		goto Eighty
 	}
@@ -1511,7 +1507,7 @@ Sixty:
 	zn = complex(ZNR, ZNI)
 	c1 = complex(C1R, C1I)
 	c2 = complex(C2R, C2I)
-	c1, c2, NW, IUF = Zs1s2(zn, c1, c2, ASCLE, ALIM, IUF)
+	c1, c2, NW, _ = Zs1s2(zn, c1, c2, ASCLE, ALIM, IUF)
 	C1R = real(c1)
 	C1I = imag(c1)
 	C2R = real(c2)
@@ -1533,8 +1529,8 @@ Eighty:
 // MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
 // REGION CABS(Z)>MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL return.
 // NZ<0 INDICATES AN OVERFLOW ON KODE=1.
-func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) (
-	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) {
+func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, RL, TOL, ELIM, ALIM float64) (
+	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, RLout, TOLout, ELIMout, ALIMout float64) {
 	var AA, AEZ, AK, AK1I, AK1R, ARG, ARM, ATOL,
 		AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
 		CZR, DFNU, DKI, DKR, DNU2, EZI, EZR, FDN, PI, P1I,
@@ -1552,7 +1548,6 @@ func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, RL, TOL,
 	CONER = 1
 	CONEI = 0
 
-	NZ = 0
 	AZ = cmplx.Abs(complex(ZR, ZI))
 	ARM = 1.0e3 * dmach[1]
 	RTR1 = math.Sqrt(ARM)
@@ -1735,8 +1730,8 @@ OneTen:
 
 // ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z)>=0.0 BY THE
 // MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
-func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL float64) (
-	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout float64) {
+func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, TOL float64) (
+	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, TOLout float64) {
 	var ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
 		CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, PTI, PTR, P1I,
 		P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
@@ -1749,7 +1744,6 @@ func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL float
 	CONEI = 0
 
 	SCLE = dmach[1] / TOL
-	NZ = 0
 	AZ = cmplx.Abs(complex(ZR, ZI))
 	IAZ = int(float32(AZ))
 	IFNU = int(float32(FNU))
diff --git a/mathext/internal/cephes/incbi.go b/mathext/internal/cephes/incbi.go
index b878af09..2b612d83 100644
--- a/mathext/internal/cephes/incbi.go
+++ b/mathext/internal/cephes/incbi.go
@@ -16,7 +16,6 @@ func Incbi(aa, bb, yy0 float64) float64 {
 	var a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt float64
 	var i, rflg, dir, nflg int
 
-	i = 0
 	if yy0 <= 0 {
 		return (0.0)
 	}
diff --git a/optimize/lbfgs.go b/optimize/lbfgs.go
index 02359dae..588a1523 100644
--- a/optimize/lbfgs.go
+++ b/optimize/lbfgs.go
@@ -70,7 +70,6 @@ func (l *LBFGS) Init(dim, tasks int) int {
 func (l *LBFGS) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	l.status, l.err = localOptimizer{}.run(l, l.GradStopThreshold, operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (l *LBFGS) initLocal(loc *Location) (Operation, error) {
diff --git a/optimize/neldermead.go b/optimize/neldermead.go
index 4615bf48..d92db09e 100644
--- a/optimize/neldermead.go
+++ b/optimize/neldermead.go
@@ -105,7 +105,6 @@ func (n *NelderMead) Init(dim, tasks int) int {
 func (n *NelderMead) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	n.status, n.err = localOptimizer{}.run(n, math.NaN(), operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (n *NelderMead) initLocal(loc *Location) (Operation, error) {
diff --git a/optimize/newton.go b/optimize/newton.go
index 9519baf3..cc8d53b1 100644
--- a/optimize/newton.go
+++ b/optimize/newton.go
@@ -82,7 +82,6 @@ func (n *Newton) Init(dim, tasks int) int {
 func (n *Newton) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	n.status, n.err = localOptimizer{}.run(n, n.GradStopThreshold, operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (n *Newton) initLocal(loc *Location) (Operation, error) {
diff --git a/stat/combin/combin.go b/stat/combin/combin.go
index 28e312b5..c349de4c 100644
--- a/stat/combin/combin.go
+++ b/stat/combin/combin.go
@@ -284,16 +284,10 @@ func IndexToCombination(dst []int, idx, n, k int) []int {
 //  [ 1 2 0 ]
 // Cartesian panics if any of the provided lengths are less than 1.
 func Cartesian(lens []int) [][]int {
-	if len(lens) == 0 {
+	rows := Card(lens)
+	if rows == 0 {
 		panic("combin: empty lengths")
 	}
-	rows := 1
-	for _, v := range lens {
-		if v < 1 {
-			panic("combin: length less than zero")
-		}
-		rows *= v
-	}
 	out := make([][]int, rows)
 	for i := 0; i < rows; i++ {
 		out[i] = SubFor(nil, i, lens)
@@ -301,6 +295,56 @@ func Cartesian(lens []int) [][]int {
 	return out
 }
 
+// Card computes the cardinality of the multi-dimensional space whose dimensions have size specified by dims
+// All length values must be positive, otherwise this will panic.
+func Card(dims []int) int {
+	if len(dims) == 0 {
+		return 0
+	}
+	card := 1
+	for _, v := range dims {
+		if v < 1 {
+			panic("combin: length less than zero")
+		}
+		card *= v
+	}
+	return card
+}
+
+// NewCartesianGenerator returns a CartesianGenerator for iterating over cartesian products which are generated on the fly.
+// All values in lens must be positive, otherwise this will panic.
+func NewCartesianGenerator(lens []int) *CartesianGenerator {
+	return &CartesianGenerator{
+		lens: lens,
+		rows: Card(lens),
+		idx:  -1,
+	}
+}
+
+// CartesianGenerator iterates over a cartesian product set.
+type CartesianGenerator struct {
+	lens []int
+	rows int
+	idx  int
+}
+
+// Next moves to the next product of the cartesian set.
+// It returns false if the generator reached the end of the cartesian set end.
+func (g *CartesianGenerator) Next() bool {
+	if g.idx+1 < g.rows {
+		g.idx++
+		return true
+	}
+	g.idx = g.rows
+	return false
+}
+
+// Product generates one product of the cartesian set according to the current index which is increased by Next().
+// Next needs to be called at least one time before this method, otherwise it will panic.
+func (g *CartesianGenerator) Product() []int {
+	return SubFor(nil, g.idx, g.lens)
+}
+
 // IdxFor converts a multi-dimensional index into a linear index for a
 // multi-dimensional space. sub specifies the index for each dimension, and dims
 // specifies the size of each dimension. IdxFor is the inverse of SubFor.
diff --git a/stat/mds/mds.go b/stat/mds/mds.go
index 5a763ea2..2eab311f 100644
--- a/stat/mds/mds.go
+++ b/stat/mds/mds.go
@@ -14,21 +14,19 @@ import (
 // TorgersonScaling converts a dissimilarity matrix to a matrix containing
 // Euclidean coordinates. TorgersonScaling places the coordinates in dst and
 // returns it and the number of positive Eigenvalues if successful.
-// If the scaling is not successful, dst is returned, but will not be a valid scaling.
-// When the scaling is successful, mds will be resized to k columns wide.
+// If the scaling is not successful, dst will be empty upon return.
+// When the scaling is successful, dst will be resized to k columns wide.
 // Eigenvalues will be copied into eigdst and returned as eig if it is provided.
 //
-// If dst is nil, a new mat.Dense is allocated. If dst is not empty,
-// the dimensions of dst and dis must match otherwise TorgersonScaling will panic.
-// The dis matrix must be square or TorgersonScaling will panic.
-func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k int, mds *mat.Dense, eig []float64) {
+// TorgersonScaling will panic if dst is not empty.
+func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k int, eig []float64) {
 	// https://doi.org/10.1007/0-387-28981-X_12
 
 	n := dis.Symmetric()
-	if dst == nil {
-		dst = mat.NewDense(n, n, nil)
-	} else if r, c := dst.Dims(); !dst.IsEmpty() && (r != n || c != n) {
-		panic(mat.ErrShape)
+	if dst.IsEmpty() {
+		dst.ReuseAs(n, n)
+	} else {
+		panic("mds: receiver matrix not empty")
 	}
 
 	b := mat.NewSymDense(n, nil)
@@ -57,7 +55,7 @@ func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k in
 	var ed mat.EigenSym
 	ok := ed.Factorize(b, true)
 	if !ok {
-		return 0, dst, eigdst
+		return 0, eigdst
 	}
 	ed.VectorsTo(dst)
 	vals := ed.Values(nil)
@@ -75,10 +73,10 @@ func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k in
 
 	var tmp mat.Dense
 	tmp.Mul(dst, mat.NewDiagonalRect(n, k, vals[:k]))
-	dst = dst.Slice(0, n, 0, k).(*mat.Dense)
+	*dst = *dst.Slice(0, n, 0, k).(*mat.Dense)
 	dst.Copy(&tmp)
 
-	return k, dst, eigdst
+	return k, eigdst
 }
 
 func reverse(values []float64, vectors blas64.General) {
diff --git a/stat/pca_cca.go b/stat/pca_cca.go
index f8715006..25158cea 100644
--- a/stat/pca_cca.go
+++ b/stat/pca_cca.go
@@ -48,24 +48,23 @@ func (c *PC) PrincipalComponents(a mat.Matrix, weights []float64) (ok bool) {
 
 // VectorsTo returns the component direction vectors of a principal components
 // analysis. The vectors are returned in the columns of a d×min(n, d) matrix.
-// If dst is not nil it must either be empty or be a d×min(n, d) matrix.
-// dst will  be used as the destination for the direction vector data. If dst
-// is nil, a new mat.Dense is allocated for the destination.
-func (c *PC) VectorsTo(dst *mat.Dense) *mat.Dense {
+//
+// If dst is empty, VectorsTo will resize dst to be d×min(n, d). When dst is
+// non-empty, VectorsTo will panic if dst is not d×min(n, d). VectorsTo will also
+// panic if the receiver does not contain a successful PC.
+func (c *PC) VectorsTo(dst *mat.Dense) {
 	if !c.ok {
 		panic("stat: use of unsuccessful principal components analysis")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); !dst.IsEmpty() && (d != c.d || n != min(c.n, c.d)) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.d, min(c.n, c.d))
+	} else {
+		if d, n := dst.Dims(); d != c.d || n != min(c.n, c.d) {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.svd.VTo(dst)
-	return dst
 }
 
 // VarsTo returns the column variances of the principal component scores,
@@ -222,26 +221,25 @@ func (c *CC) CorrsTo(dst []float64) []float64 {
 // LeftTo returns the left eigenvectors of the canonical correlation matrix if
 // spheredSpace is true. If spheredSpace is false it returns these eigenvectors
 // back-transformed to the original data space.
-// If dst is not nil it must either be empty or be an xd×yd matrix where xd
-// and yd are the number of variables in the input x and y matrices. dst will
-// be used as the destination for the vector data. If dst is nil, a new
-// mat.Dense is allocated for the destination.
-func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
+//
+// If dst is empty, LeftTo will resize dst to be xd×yd. When dst is
+// non-empty, LeftTo will panic if dst is not xd×yd. LeftTo will also
+// panic if the receiver does not contain a successful CC.
+func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) {
 	if !c.ok || c.n < 2 {
 		panic("stat: canonical correlations missing or invalid")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); !dst.IsEmpty() && (n != c.yd || d != c.xd) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.xd, c.yd)
+	} else {
+		if d, n := dst.Dims(); d != c.xd || n != c.yd {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.c.UTo(dst)
 	if spheredSpace {
-		return dst
+		return
 	}
 
 	xs := c.x.Values(nil)
@@ -252,32 +250,30 @@ func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
 
 	dst.Product(xv, xv.T(), dst)
 	dst.Scale(math.Sqrt(float64(c.n-1)), dst)
-	return dst
 }
 
 // RightTo returns the right eigenvectors of the canonical correlation matrix if
 // spheredSpace is true. If spheredSpace is false it returns these eigenvectors
 // back-transformed to the original data space.
-// If dst is not nil it must either be zero-sized or be an yd×yd matrix where yd
-// is the number of variables in the input y matrix. dst will
-// be used as the destination for the vector data. If dst is nil, a new
-// mat.Dense is allocated for the destination.
-func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
+//
+// If dst is empty, RightTo will resize dst to be yd×yd. When dst is
+// non-empty, RightTo will panic if dst is not yd×yd. RightTo will also
+// panic if the receiver does not contain a successful CC.
+func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) {
 	if !c.ok || c.n < 2 {
 		panic("stat: canonical correlations missing or invalid")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); (n != 0 || d != 0) && (n != c.yd || d != c.yd) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.yd, c.yd)
+	} else {
+		if d, n := dst.Dims(); d != c.yd || n != c.yd {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.c.VTo(dst)
 	if spheredSpace {
-		return dst
+		return
 	}
 
 	ys := c.y.Values(nil)
@@ -288,7 +284,6 @@ func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
 
 	dst.Product(yv, yv.T(), dst)
 	dst.Scale(math.Sqrt(float64(c.n-1)), dst)
-	return dst
 }
 
 func svdFactorizeCentered(work *mat.SVD, m mat.Matrix, weights []float64) (svd *mat.SVD, ok bool) {


Test patch:
diff --git a/graph/implicit_example_test.go b/graph/implicit_example_test.go
new file mode 100644
index 00000000..25c3e4d0
--- /dev/null
+++ b/graph/implicit_example_test.go
@@ -0,0 +1,385 @@
+// Copyright ©2014 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 graph_test
+
+import (
+	"fmt"
+
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/iterator"
+	"gonum.org/v1/gonum/graph/simple"
+	"gonum.org/v1/gonum/graph/topo"
+)
+
+// GraphNode is a node in an implicit graph.
+type GraphNode struct {
+	id        int64
+	neighbors []graph.Node
+	roots     []*GraphNode
+}
+
+// NewGraphNode returns a new GraphNode.
+func NewGraphNode(id int64) *GraphNode {
+	return &GraphNode{id: id}
+}
+
+// Node allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) Node(id int64) graph.Node {
+	if id == g.id {
+		return g
+	}
+
+	seen := map[int64]struct{}{g.id: {}}
+	for _, root := range g.roots {
+		if root.ID() == id {
+			return root
+		}
+
+		if root.has(seen, id) {
+			return root
+		}
+	}
+
+	for _, n := range g.neighbors {
+		if n.ID() == id {
+			return n
+		}
+
+		if gn, ok := n.(*GraphNode); ok {
+			if gn.has(seen, id) {
+				return gn
+			}
+		}
+	}
+
+	return nil
+}
+
+func (g *GraphNode) has(seen map[int64]struct{}, id int64) bool {
+	for _, root := range g.roots {
+		if _, ok := seen[root.ID()]; ok {
+			continue
+		}
+
+		seen[root.ID()] = struct{}{}
+		if root.ID() == id {
+			return true
+		}
+
+		if root.has(seen, id) {
+			return true
+		}
+
+	}
+
+	for _, n := range g.neighbors {
+		if _, ok := seen[n.ID()]; ok {
+			continue
+		}
+
+		seen[n.ID()] = struct{}{}
+		if n.ID() == id {
+			return true
+		}
+
+		if gn, ok := n.(*GraphNode); ok {
+			if gn.has(seen, id) {
+				return true
+			}
+		}
+	}
+
+	return false
+}
+
+// Nodes allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) Nodes() graph.Nodes {
+	nodes := []graph.Node{g}
+	seen := map[int64]struct{}{g.id: {}}
+
+	for _, root := range g.roots {
+		nodes = append(nodes, root)
+		seen[root.ID()] = struct{}{}
+
+		nodes = root.nodes(nodes, seen)
+	}
+
+	for _, n := range g.neighbors {
+		nodes = append(nodes, n)
+		seen[n.ID()] = struct{}{}
+
+		if gn, ok := n.(*GraphNode); ok {
+			nodes = gn.nodes(nodes, seen)
+		}
+	}
+
+	return iterator.NewOrderedNodes(nodes)
+}
+
+func (g *GraphNode) nodes(dst []graph.Node, seen map[int64]struct{}) []graph.Node {
+	for _, root := range g.roots {
+		if _, ok := seen[root.ID()]; ok {
+			continue
+		}
+		seen[root.ID()] = struct{}{}
+		dst = append(dst, graph.Node(root))
+
+		dst = root.nodes(dst, seen)
+	}
+
+	for _, n := range g.neighbors {
+		if _, ok := seen[n.ID()]; ok {
+			continue
+		}
+
+		dst = append(dst, n)
+		if gn, ok := n.(*GraphNode); ok {
+			dst = gn.nodes(dst, seen)
+		}
+	}
+
+	return dst
+}
+
+// From allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) From(id int64) graph.Nodes {
+	if id == g.ID() {
+		return iterator.NewOrderedNodes(g.neighbors)
+	}
+
+	seen := map[int64]struct{}{g.id: {}}
+	for _, root := range g.roots {
+		seen[root.ID()] = struct{}{}
+
+		if result := root.findNeighbors(id, seen); result != nil {
+			return iterator.NewOrderedNodes(result)
+		}
+	}
+
+	for _, n := range g.neighbors {
+		seen[n.ID()] = struct{}{}
+
+		if gn, ok := n.(*GraphNode); ok {
+			if result := gn.findNeighbors(id, seen); result != nil {
+				return iterator.NewOrderedNodes(result)
+			}
+		}
+	}
+
+	return nil
+}
+
+func (g *GraphNode) findNeighbors(id int64, seen map[int64]struct{}) []graph.Node {
+	if id == g.ID() {
+		return g.neighbors
+	}
+
+	for _, root := range g.roots {
+		if _, ok := seen[root.ID()]; ok {
+			continue
+		}
+		seen[root.ID()] = struct{}{}
+
+		if result := root.findNeighbors(id, seen); result != nil {
+			return result
+		}
+	}
+
+	for _, n := range g.neighbors {
+		if _, ok := seen[n.ID()]; ok {
+			continue
+		}
+		seen[n.ID()] = struct{}{}
+
+		if gn, ok := n.(*GraphNode); ok {
+			if result := gn.findNeighbors(id, seen); result != nil {
+				return result
+			}
+		}
+	}
+
+	return nil
+}
+
+// HasEdgeBetween allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) HasEdgeBetween(uid, vid int64) bool {
+	return g.EdgeBetween(uid, vid) != nil
+}
+
+// Edge allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) Edge(uid, vid int64) graph.Edge {
+	return g.EdgeBetween(uid, vid)
+}
+
+// EdgeBetween allows GraphNode to satisfy the graph.Graph interface.
+func (g *GraphNode) EdgeBetween(uid, vid int64) graph.Edge {
+	if uid == g.id || vid == g.id {
+		for _, n := range g.neighbors {
+			if n.ID() == uid || n.ID() == vid {
+				return simple.Edge{F: g, T: n}
+			}
+		}
+		return nil
+	}
+
+	seen := map[int64]struct{}{g.id: {}}
+	for _, root := range g.roots {
+		seen[root.ID()] = struct{}{}
+		if result := root.edgeBetween(uid, vid, seen); result != nil {
+			return result
+		}
+	}
+
+	for _, n := range g.neighbors {
+		seen[n.ID()] = struct{}{}
+		if gn, ok := n.(*GraphNode); ok {
+			if result := gn.edgeBetween(uid, vid, seen); result != nil {
+				return result
+			}
+		}
+	}
+
+	return nil
+}
+
+func (g *GraphNode) edgeBetween(uid, vid int64, seen map[int64]struct{}) graph.Edge {
+	if uid == g.id || vid == g.id {
+		for _, n := range g.neighbors {
+			if n.ID() == uid || n.ID() == vid {
+				return simple.Edge{F: g, T: n}
+			}
+		}
+		return nil
+	}
+
+	for _, root := range g.roots {
+		if _, ok := seen[root.ID()]; ok {
+			continue
+		}
+		seen[root.ID()] = struct{}{}
+		if result := root.edgeBetween(uid, vid, seen); result != nil {
+			return result
+		}
+	}
+
+	for _, n := range g.neighbors {
+		if _, ok := seen[n.ID()]; ok {
+			continue
+		}
+
+		seen[n.ID()] = struct{}{}
+		if gn, ok := n.(*GraphNode); ok {
+			if result := gn.edgeBetween(uid, vid, seen); result != nil {
+				return result
+			}
+		}
+	}
+
+	return nil
+}
+
+// ID allows GraphNode to satisfy the graph.Node interface.
+func (g *GraphNode) ID() int64 {
+	return g.id
+}
+
+// AddMeighbor adds an edge between g and n.
+func (g *GraphNode) AddNeighbor(n *GraphNode) {
+	g.neighbors = append(g.neighbors, graph.Node(n))
+}
+
+// AddRoot adds provides an entrance into the graph g from n.
+func (g *GraphNode) AddRoot(n *GraphNode) {
+	g.roots = append(g.roots, n)
+}
+
+func Example_implicit() {
+	// This example shows the construction of the following graph
+	// using the implicit graph type above.
+	//
+	// The visual representation of the graph can be seen at
+	// https://graphviz.gitlab.io/_pages/Gallery/undirected/fdpclust.html
+	//
+	// graph G {
+	// 	e
+	// 	subgraph clusterA {
+	// 		a -- b
+	// 		subgraph clusterC {
+	// 			C -- D
+	// 		}
+	// 	}
+	// 	subgraph clusterB {
+	// 		d -- f
+	// 	}
+	// 	d -- D
+	// 	e -- clusterB
+	// 	clusterC -- clusterB
+	// }
+
+	// graph G {
+	G := NewGraphNode(0)
+	// 	e
+	e := NewGraphNode(1)
+
+	// 	subgraph clusterA {
+	clusterA := NewGraphNode(2)
+	// 		a -- b
+	a := NewGraphNode(3)
+	b := NewGraphNode(4)
+	a.AddNeighbor(b)
+	b.AddNeighbor(a)
+	clusterA.AddRoot(a)
+	clusterA.AddRoot(b)
+
+	// 		subgraph clusterC {
+	clusterC := NewGraphNode(5)
+	// 			C -- D
+	C := NewGraphNode(6)
+	D := NewGraphNode(7)
+	C.AddNeighbor(D)
+	D.AddNeighbor(C)
+
+	clusterC.AddRoot(C)
+	clusterC.AddRoot(D)
+	// 		}
+	clusterA.AddRoot(clusterC)
+	// 	}
+
+	// 	subgraph clusterB {
+	clusterB := NewGraphNode(8)
+	// 		d -- f
+	d := NewGraphNode(9)
+	f := NewGraphNode(10)
+	d.AddNeighbor(f)
+	f.AddNeighbor(d)
+	clusterB.AddRoot(d)
+	clusterB.AddRoot(f)
+	// 	}
+
+	// 	d -- D
+	d.AddNeighbor(D)
+	D.AddNeighbor(d)
+
+	// 	e -- clusterB
+	e.AddNeighbor(clusterB)
+	clusterB.AddNeighbor(e)
+
+	// 	clusterC -- clusterB
+	clusterC.AddNeighbor(clusterB)
+	clusterB.AddNeighbor(clusterC)
+
+	G.AddRoot(e)
+	G.AddRoot(clusterA)
+	G.AddRoot(clusterB)
+	// }
+
+	if topo.IsPathIn(G, []graph.Node{C, D, d, f}) {
+		fmt.Println("C--D--d--f is a path in G.")
+	}
+
+	// Output:
+	//
+	// C--D--d--f is a path in G.
+}
diff --git a/graph/iterator/nodes_test.go b/graph/iterator/nodes_test.go
index f7dab387..d8db2ada 100644
--- a/graph/iterator/nodes_test.go
+++ b/graph/iterator/nodes_test.go
@@ -97,3 +97,49 @@ func TestImplicitNodesIterate(t *testing.T) {
 		}
 	}
 }
+
+var nodesTests = []struct {
+	nodes map[int64]graph.Node
+}{
+	{nodes: nil},
+	{nodes: make(map[int64]graph.Node)},
+	{nodes: map[int64]graph.Node{1: simple.Node(1)}},
+	{nodes: map[int64]graph.Node{1: simple.Node(1), 2: simple.Node(2), 3: simple.Node(3), 5: simple.Node(5)}},
+	{nodes: map[int64]graph.Node{5: simple.Node(5), 3: simple.Node(3), 2: simple.Node(2), 1: simple.Node(1)}},
+}
+
+func TestNodesIterate(t *testing.T) {
+	for _, test := range nodesTests {
+		it := iterator.NewNodes(test.nodes)
+		for i := 0; i < 2; i++ {
+			if it.Len() != len(test.nodes) {
+				t.Errorf("unexpected iterator length for round %d: got:%d want:%d", i, it.Len(), len(test.nodes))
+			}
+			var got map[int64]graph.Node
+			if test.nodes != nil {
+				got = make(map[int64]graph.Node)
+			}
+			for it.Next() {
+				n := it.Node()
+				got[n.ID()] = n
+				if len(got)+it.Len() != len(test.nodes) {
+					t.Errorf("unexpected iterator length during iteration for round %d: got:%d want:%d", i, it.Len(), len(test.nodes))
+				}
+			}
+			want := test.nodes
+			if !reflect.DeepEqual(got, want) {
+				t.Errorf("unexpected iterator output for round %d: got:%#v want:%#v", i, got, want)
+			}
+			func() {
+				defer func() {
+					r := recover()
+					if r != nil {
+						t.Errorf("unexpected panic: %v", r)
+					}
+				}()
+				it.Next()
+			}()
+			it.Reset()
+		}
+	}
+}
diff --git a/graph/iterator/nodes_unordered_test.go b/graph/iterator/nodes_unordered_test.go
deleted file mode 100644
index bace9684..00000000
--- a/graph/iterator/nodes_unordered_test.go
+++ /dev/null
@@ -1,64 +0,0 @@
-// Copyright ©2018 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.
-
-// +build go1.12
-
-// TODO(kortschak): Merge this into nodes_test.go when go1.11 is no longer supported.
-
-package iterator_test
-
-import (
-	"reflect"
-	"testing"
-
-	"gonum.org/v1/gonum/graph"
-	"gonum.org/v1/gonum/graph/iterator"
-	"gonum.org/v1/gonum/graph/simple"
-)
-
-var nodesTests = []struct {
-	nodes map[int64]graph.Node
-}{
-	{nodes: nil},
-	{nodes: make(map[int64]graph.Node)},
-	{nodes: map[int64]graph.Node{1: simple.Node(1)}},
-	{nodes: map[int64]graph.Node{1: simple.Node(1), 2: simple.Node(2), 3: simple.Node(3), 5: simple.Node(5)}},
-	{nodes: map[int64]graph.Node{5: simple.Node(5), 3: simple.Node(3), 2: simple.Node(2), 1: simple.Node(1)}},
-}
-
-func TestNodesIterate(t *testing.T) {
-	for _, test := range nodesTests {
-		it := iterator.NewNodes(test.nodes)
-		for i := 0; i < 2; i++ {
-			if it.Len() != len(test.nodes) {
-				t.Errorf("unexpected iterator length for round %d: got:%d want:%d", i, it.Len(), len(test.nodes))
-			}
-			var got map[int64]graph.Node
-			if test.nodes != nil {
-				got = make(map[int64]graph.Node)
-			}
-			for it.Next() {
-				n := it.Node()
-				got[n.ID()] = n
-				if len(got)+it.Len() != len(test.nodes) {
-					t.Errorf("unexpected iterator length during iteration for round %d: got:%d want:%d", i, it.Len(), len(test.nodes))
-				}
-			}
-			want := test.nodes
-			if !reflect.DeepEqual(got, want) {
-				t.Errorf("unexpected iterator output for round %d: got:%#v want:%#v", i, got, want)
-			}
-			func() {
-				defer func() {
-					r := recover()
-					if r != nil {
-						t.Errorf("unexpected panic: %v", r)
-					}
-				}()
-				it.Next()
-			}()
-			it.Reset()
-		}
-	}
-}
diff --git a/graph/path/dynamic/dumper_test.go b/graph/path/dynamic/dumper_test.go
index 1ab77210..05646954 100644
--- a/graph/path/dynamic/dumper_test.go
+++ b/graph/path/dynamic/dumper_test.go
@@ -96,7 +96,7 @@ func (d *dumper) dump(withpath bool) {
 					if n.g != n.rhs {
 						fmt.Fprintf(w, "key:%.3f", n.key)
 					}
-					if n.key == n.key {
+					if !n.key.isBadKey() {
 						// Mark keys for nodes in the priority queue.
 						// We use NaN inequality for this check since all
 						// keys not in the queue must have their key set
diff --git a/integrate/romberg_test.go b/integrate/romberg_test.go
new file mode 100644
index 00000000..00c34adc
--- /dev/null
+++ b/integrate/romberg_test.go
@@ -0,0 +1,65 @@
+// Copyright ©2019 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 integrate
+
+import (
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/integrate/testquad"
+)
+
+func TestRomberg(t *testing.T) {
+	for i, test := range []struct {
+		integral testquad.Integral
+		n        int
+		tol      float64
+	}{
+		{integral: testquad.Constant(0), n: 3, tol: 0},
+		{integral: testquad.Constant(0), n: 1<<5 + 1, tol: 0},
+		{integral: testquad.Poly(0), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(0), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Poly(2), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(2), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Poly(3), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(3), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Poly(4), n: 5, tol: 1e-14},
+		{integral: testquad.Poly(4), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Poly(5), n: 5, tol: 1e-14},
+		{integral: testquad.Poly(5), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Sin(), n: 1<<3 + 1, tol: 1e-10},
+		{integral: testquad.Sin(), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.XExpMinusX(), n: 1<<3 + 1, tol: 1e-9},
+		{integral: testquad.XExpMinusX(), n: 1<<5 + 1, tol: 1e-14},
+		{integral: testquad.Sqrt(), n: 1<<10 + 1, tol: 1e-5},
+		{integral: testquad.ExpOverX2Plus1(), n: 1<<4 + 1, tol: 1e-7},
+		{integral: testquad.ExpOverX2Plus1(), n: 1<<6 + 1, tol: 1e-14},
+	} {
+		n := test.n
+		a := test.integral.A
+		b := test.integral.B
+
+		x := make([]float64, n)
+		floats.Span(x, a, b)
+
+		y := make([]float64, n)
+		for i, xi := range x {
+			y[i] = test.integral.F(xi)
+		}
+
+		dx := (b - a) / float64(n-1)
+		got := Romberg(y, dx)
+
+		want := test.integral.Value
+		diff := math.Abs(got - want)
+		if diff > test.tol {
+			t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+				i, test.integral.Name, n, got, want, diff)
+		}
+	}
+}
diff --git a/integrate/simpsons_test.go b/integrate/simpsons_test.go
new file mode 100644
index 00000000..cb6dc543
--- /dev/null
+++ b/integrate/simpsons_test.go
@@ -0,0 +1,58 @@
+// Copyright ©2019 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 integrate
+
+import (
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/integrate/testquad"
+)
+
+func TestSimpsons(t *testing.T) {
+	rnd := rand.New(rand.NewSource(1))
+	for i, test := range []struct {
+		integral testquad.Integral
+		n        int
+		tol      float64
+	}{
+		{integral: testquad.Constant(0), n: 3, tol: 0},
+		{integral: testquad.Constant(0), n: 10, tol: 0},
+		{integral: testquad.Poly(0), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(0), n: 10, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 10, tol: 1e-14},
+		{integral: testquad.Poly(2), n: 3, tol: 1e-14},
+		{integral: testquad.Poly(2), n: 10, tol: 1e-14},
+		{integral: testquad.Poly(3), n: 1e3, tol: 1e-8},
+		{integral: testquad.Poly(4), n: 1e3, tol: 1e-8},
+		{integral: testquad.Poly(5), n: 1e3, tol: 1e-7},
+		{integral: testquad.Sin(), n: 1e2, tol: 1e-8},
+		{integral: testquad.XExpMinusX(), n: 1e2, tol: 1e-8},
+		{integral: testquad.Sqrt(), n: 1e4, tol: 1e-6},
+		{integral: testquad.ExpOverX2Plus1(), n: 1e2, tol: 1e-7},
+	} {
+		n := test.n
+		a := test.integral.A
+		b := test.integral.B
+
+		x := jitterSpan(n, a, b, rnd)
+		y := make([]float64, n)
+		for i, xi := range x {
+			y[i] = test.integral.F(xi)
+		}
+
+		got := Simpsons(x, y)
+
+		want := test.integral.Value
+		diff := math.Abs(got - want)
+		if diff > test.tol {
+			t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+				i, test.integral.Name, n, got, want, diff)
+		}
+	}
+}
diff --git a/integrate/testquad/testquad.go b/integrate/testquad/testquad.go
new file mode 100644
index 00000000..f6b6a8f6
--- /dev/null
+++ b/integrate/testquad/testquad.go
@@ -0,0 +1,106 @@
+// Copyright ©2019 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 testquad provides integrals for testing quadrature algorithms.
+package testquad
+
+import (
+	"fmt"
+	"math"
+)
+
+// Integral is a definite integral
+//  ∫_a^b f(x)dx
+// with a known value.
+type Integral struct {
+	Name  string
+	A, B  float64               // Integration limits
+	F     func(float64) float64 // Integrand
+	Value float64
+}
+
+// Constant returns the integral of a constant function
+//  ∫_{-1}^2 alpha dx
+func Constant(alpha float64) Integral {
+	return Integral{
+		Name: fmt.Sprintf("∫_{-1}^{2} %vdx", alpha),
+		A:    -1,
+		B:    2,
+		F: func(float64) float64 {
+			return alpha
+		},
+		Value: 3 * alpha,
+	}
+}
+
+// Poly returns the integral of a polynomial
+//  ∫_{-1}^2 x^degree dx
+func Poly(degree int) Integral {
+	d := float64(degree)
+	return Integral{
+		Name: fmt.Sprintf("∫_{-1}^{2} x^%vdx", degree),
+		A:    -1,
+		B:    2,
+		F: func(x float64) float64 {
+			return math.Pow(x, d)
+		},
+		Value: (math.Pow(2, d+1) - math.Pow(-1, d+1)) / (d + 1),
+	}
+}
+
+// Sin returns the integral
+//  ∫_0^1 sin(x)dx
+func Sin() Integral {
+	return Integral{
+		Name: "∫_0^1 sin(x)dx",
+		A:    0,
+		B:    1,
+		F: func(x float64) float64 {
+			return math.Sin(x)
+		},
+		Value: 1 - math.Cos(1),
+	}
+}
+
+// XExpMinusX returns the integral
+//  ∫_0^1 x*exp(-x)dx
+func XExpMinusX() Integral {
+	return Integral{
+		Name: "∫_0^1 x*exp(-x)dx",
+		A:    0,
+		B:    1,
+		F: func(x float64) float64 {
+			return x * math.Exp(-x)
+		},
+		Value: (math.E - 2) / math.E,
+	}
+}
+
+// Sqrt returns the integral
+//  ∫_0^1 sqrt(x)dx
+func Sqrt() Integral {
+	return Integral{
+		Name: "∫_0^1 sqrt(x)dx",
+		A:    0,
+		B:    1,
+		F: func(x float64) float64 {
+			return math.Sqrt(x)
+		},
+		Value: 2 / 3.0,
+	}
+}
+
+// ExpOverX2Plus1 returns the integral
+//  ∫_0^1 exp(x)/(x*x+1)dx
+func ExpOverX2Plus1() Integral {
+	return Integral{
+		Name: "∫_0^1 exp(x)/(x*x+1)dx",
+		A:    0,
+		B:    1,
+		F: func(x float64) float64 {
+			return math.Exp(x) / (x*x + 1)
+		},
+		Value: 1.270724139833620220138,
+	}
+}
diff --git a/integrate/trapezoidal_test.go b/integrate/trapezoidal_test.go
index 1a4451ee..68b73792 100644
--- a/integrate/trapezoidal_test.go
+++ b/integrate/trapezoidal_test.go
@@ -8,75 +8,66 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/integrate/testquad"
 )
 
 func TestTrapezoidal(t *testing.T) {
-	const N = 1e6
-	x := floats.Span(make([]float64, N), 0, 1)
+	rnd := rand.New(rand.NewSource(1))
 	for i, test := range []struct {
-		x    []float64
-		f    func(x float64) float64
-		want float64
+		integral testquad.Integral
+		n        int
+		tol      float64
 	}{
-		{
-			x:    x,
-			f:    func(x float64) float64 { return x },
-			want: 0.5,
-		},
-		{
-			x:    floats.Span(make([]float64, N), -1, 1),
-			f:    func(x float64) float64 { return x },
-			want: 0,
-		},
-		{
-			x:    x,
-			f:    func(x float64) float64 { return x + 10 },
-			want: 10.5,
-		},
-		{
-			x:    x,
-			f:    func(x float64) float64 { return 3*x*x + 10 },
-			want: 11,
-		},
-		{
-			x:    x,
-			f:    func(x float64) float64 { return math.Exp(x) },
-			want: 1.7182818284591876,
-		},
-		{
-			x:    floats.Span(make([]float64, N), 0, math.Pi),
-			f:    func(x float64) float64 { return math.Cos(x) },
-			want: 0,
-		},
-		{
-			x:    floats.Span(make([]float64, N), 0, 2*math.Pi),
-			f:    func(x float64) float64 { return math.Cos(x) },
-			want: 0,
-		},
-		{
-			x:    floats.Span(make([]float64, N*10), 0, math.Pi),
-			f:    func(x float64) float64 { return math.Sin(x) },
-			want: 2,
-		},
-		{
-			x:    floats.Span(make([]float64, N*10), 0, 0.5*math.Pi),
-			f:    func(x float64) float64 { return math.Sin(x) },
-			want: 1,
-		},
-		{
-			x:    floats.Span(make([]float64, N), 0, 2*math.Pi),
-			f:    func(x float64) float64 { return math.Sin(x) },
-			want: 0,
-		},
+		{integral: testquad.Constant(0), n: 2, tol: 0},
+		{integral: testquad.Constant(0), n: 10, tol: 0},
+		{integral: testquad.Poly(0), n: 2, tol: 1e-14},
+		{integral: testquad.Poly(0), n: 10, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 2, tol: 1e-14},
+		{integral: testquad.Poly(1), n: 10, tol: 1e-14},
+		{integral: testquad.Poly(2), n: 1e5, tol: 1e-8},
+		{integral: testquad.Poly(3), n: 1e5, tol: 1e-8},
+		{integral: testquad.Poly(4), n: 1e5, tol: 1e-7},
+		{integral: testquad.Poly(5), n: 1e5, tol: 1e-7},
+		{integral: testquad.Sin(), n: 1e5, tol: 1e-11},
+		{integral: testquad.XExpMinusX(), n: 1e5, tol: 1e-10},
+		{integral: testquad.Sqrt(), n: 1e5, tol: 1e-8},
+		{integral: testquad.ExpOverX2Plus1(), n: 1e5, tol: 1e-10},
 	} {
-		y := make([]float64, len(test.x))
-		for i, v := range test.x {
-			y[i] = test.f(v)
+		n := test.n
+		a := test.integral.A
+		b := test.integral.B
+
+		x := jitterSpan(n, a, b, rnd)
+		y := make([]float64, n)
+		for i, xi := range x {
+			y[i] = test.integral.F(xi)
 		}
-		v := Trapezoidal(test.x, y)
-		if !floats.EqualWithinAbs(v, test.want, 1e-12) {
-			t.Errorf("test #%d: got=%v want=%v\n", i, v, test.want)
+
+		got := Trapezoidal(x, y)
+
+		want := test.integral.Value
+		diff := math.Abs(got - want)
+		if diff > test.tol {
+			t.Errorf("Test #%d: %v, n=%v: unexpected result; got=%v want=%v diff=%v",
+				i, test.integral.Name, n, got, want, diff)
 		}
 	}
 }
+
+func jitterSpan(n int, a, b float64, rnd *rand.Rand) []float64 {
+	dx := (b - a) / float64(n-1)
+	x := make([]float64, n)
+	x[0] = a
+	for i := 1; i < n-1; i++ {
+		// Set x[i] to its regular location.
+		x[i] = a + float64(i)*dx
+		// Generate a random number in [-1,1).
+		jitter := 2*rnd.Float64() - 1
+		// Jitter x[i] without crossing over its neighbors.
+		x[i] += 0.4 * jitter * dx
+	}
+	x[n-1] = b
+	return x
+}
diff --git a/mat/dense_example_test.go b/mat/dense_example_test.go
index de4a18b2..d73e1398 100644
--- a/mat/dense_example_test.go
+++ b/mat/dense_example_test.go
@@ -6,6 +6,7 @@ package mat_test
 
 import (
 	"fmt"
+	"log"
 
 	"gonum.org/v1/gonum/mat"
 )
@@ -118,7 +119,10 @@ func ExampleDense_Inverse() {
 	var ia mat.Dense
 
 	// Take the inverse of a and place the result in ia.
-	ia.Inverse(a)
+	err := ia.Inverse(a)
+	if err != nil {
+		log.Fatalf("a is not invertible: %v", err)
+	}
 
 	// Print the result using the formatter.
 	fa := mat.Formatted(&ia, mat.Prefix("     "), mat.Squeeze())
@@ -141,7 +145,10 @@ func ExampleDense_Inverse() {
 		0, 2,
 	})
 	var x mat.Dense
-	x.Solve(a, b)
+	err = x.Solve(a, b)
+	if err != nil {
+		log.Fatalf("no solution: %v", err)
+	}
 
 	// Print the result using the formatter.
 	fx := mat.Formatted(&x, mat.Prefix("    "), mat.Squeeze())
diff --git a/mathext/internal/amos/amos_test.go b/mathext/internal/amos/amos_test.go
index 3f65cdc9..fdcbbfc1 100644
--- a/mathext/internal/amos/amos_test.go
+++ b/mathext/internal/amos/amos_test.go
@@ -229,7 +229,7 @@ func zkscltest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRRamos, ZRIamos, FNUamos, Namos, YRamos, YIamos, NZamos, RZRamos, RZIamos, ASCLEamos, TOLamos, ELIMamos :=
-		Zkscl(ZRR, ZRI, FNU, n, yramos, yiamos, NZ, RZR, RZI, ASCLE, tol, ELIM)
+		Zkscl(ZRR, ZRI, FNU, n, yramos, yiamos, RZR, RZI, ASCLE, tol, ELIM)
 
 	sameF64(t, "zkscl zrr", ZRRfort, ZRRamos)
 	sameF64(t, "zkscl zri", ZRIfort, ZRIamos)
@@ -265,7 +265,7 @@ func zmlritest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, TOLamos :=
-		Zmlri(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, tol)
+		Zmlri(ZR, ZI, FNU, KODE, n, yramos, yiamos, tol)
 
 	sameF64(t, "zmlri zr", ZRfort, ZRamos)
 	sameF64(t, "zmlri zi", ZIfort, ZIamos)
@@ -358,7 +358,7 @@ func zasyitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, RLamos, TOLamos, ELIMamos, ALIMamos :=
-		Zasyi(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, RL, tol, ELIM, ALIM)
+		Zasyi(ZR, ZI, FNU, KODE, n, yramos, yiamos, RL, tol, ELIM, ALIM)
 
 	sameF64(t, "zasyi zr", ZRfort, ZRamos)
 	sameF64(t, "zasyi zr", ZIfort, ZIamos)
@@ -396,7 +396,7 @@ func zbknutest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, TOLamos, ELIMamos, ALIMamos :=
-		Zbknu(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, tol, ELIM, ALIM)
+		Zbknu(ZR, ZI, FNU, KODE, n, yramos, yiamos, tol, ELIM, ALIM)
 
 	sameF64(t, "zbknu zr", ZRfort, ZRamos)
 	sameF64(t, "zbknu zr", ZIfort, ZIamos)
@@ -418,12 +418,13 @@ func zairytest(t *testing.T, x []float64, kode, id int) {
 	KODE := kode
 	ID := id
 
-	AIRfort, AIIfort, NZfort := zairyOrig(ZR, ZI, ID, KODE)
-	AIRamos, AIIamos, NZamos := Zairy(ZR, ZI, ID, KODE)
+	AIRfort, AIIfort, NZfort, IERRfort := zairyOrig(ZR, ZI, ID, KODE)
+	AIRamos, AIIamos, NZamos, IERRamos := Zairy(ZR, ZI, ID, KODE)
 
 	sameF64Approx(t, "zairy air", AIRfort, AIRamos, 1e-12)
 	sameF64Approx(t, "zairy aii", AIIfort, AIIamos, 1e-12)
 	sameInt(t, "zairy nz", NZfort, NZamos)
+	sameInt(t, "zairy ierr", IERRfort, IERRamos)
 }
 
 func zacaitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi []float64, kode int) {
@@ -449,7 +450,7 @@ func zacaitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, MRamos, Namos, YRamos, YIamos, NZamos, RLamos, TOLamos, ELIMamos, ALIMamos :=
-		Zacai(ZR, ZI, FNU, KODE, MR, n, yramos, yiamos, NZ, RL, tol, ELIM, ALIM)
+		Zacai(ZR, ZI, FNU, KODE, MR, n, yramos, yiamos, RL, tol, ELIM, ALIM)
 
 	sameF64(t, "zacai zr", ZRfort, ZRamos)
 	sameF64(t, "zacai zi", ZIfort, ZIamos)
diff --git a/mathext/internal/amos/origcode_test.go b/mathext/internal/amos/origcode_test.go
index 56c9e798..e9cf23a8 100644
--- a/mathext/internal/amos/origcode_test.go
+++ b/mathext/internal/amos/origcode_test.go
@@ -2,6 +2,7 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
+//nolint
 package amos
 
 import (
@@ -34,7 +35,7 @@ func max0(a, b int) int {
 	return b
 }
 
-func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
+func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ, IERR int) {
 	// zairy is adapted from the original Netlib code by Donald Amos.
 	// http://www.netlib.no/netlib/amos/zairy.f
 
@@ -171,7 +172,7 @@ func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
 		DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
 		S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
 		ZEROR, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB float64
-	var IERR, IFLAG, K, K1, K2, MR, NN int
+	var IFLAG, K, K1, K2, MR, NN int
 	var tmp complex128
 
 	// Extra element for padding.
diff --git a/spatial/kdtree/kdtree_test.go b/spatial/kdtree/kdtree_test.go
index f674edd3..48cd9563 100644
--- a/spatial/kdtree/kdtree_test.go
+++ b/spatial/kdtree/kdtree_test.go
@@ -660,9 +660,9 @@ func dotFile(t *Tree, label, dotString string) (err error) {
 	}
 	defer f.Close()
 	if dotString == "" {
-		fmt.Fprintf(f, dot(t, label))
+		fmt.Fprint(f, dot(t, label))
 	} else {
-		fmt.Fprintf(f, dotString)
+		fmt.Fprint(f, dotString)
 	}
 	return
 }
diff --git a/spatial/vptree/vptree_test.go b/spatial/vptree/vptree_test.go
index 7dc81903..66e00dba 100644
--- a/spatial/vptree/vptree_test.go
+++ b/spatial/vptree/vptree_test.go
@@ -549,9 +549,9 @@ func dotFile(t *Tree, label, dotString string) (err error) {
 	}
 	defer f.Close()
 	if dotString == "" {
-		fmt.Fprintf(f, dot(t, label))
+		fmt.Fprint(f, dot(t, label))
 	} else {
-		fmt.Fprintf(f, dotString)
+		fmt.Fprint(f, dotString)
 	}
 	return
 }
diff --git a/stat/cca_example_test.go b/stat/cca_example_test.go
index 44ef9cba..a8f738e5 100644
--- a/stat/cca_example_test.go
+++ b/stat/cca_example_test.go
@@ -85,16 +85,17 @@ func ExampleCC() {
 	}
 
 	// Unpack cc.
+	var pVecs, qVecs, phiVs, psiVs mat.Dense
 	ccors := cc.CorrsTo(nil)
-	pVecs := cc.LeftTo(nil, true)
-	qVecs := cc.RightTo(nil, true)
-	phiVs := cc.LeftTo(nil, false)
-	psiVs := cc.RightTo(nil, false)
+	cc.LeftTo(&pVecs, true)
+	cc.RightTo(&qVecs, true)
+	cc.LeftTo(&phiVs, false)
+	cc.RightTo(&psiVs, false)
 
 	// Canonical Correlation Matrix, or the correlations between the sphered
 	// data.
 	var corSph mat.Dense
-	corSph.CloneFrom(pVecs)
+	corSph.CloneFrom(&pVecs)
 	col := make([]float64, xd)
 	for j := 0; j < yd; j++ {
 		mat.Col(col, j, &corSph)
@@ -110,13 +111,13 @@ func ExampleCC() {
 	fmt.Printf("\n\nccors = %.4f", ccors)
 
 	// Left and right eigenvectors of the canonical correlation matrix.
-	fmt.Printf("\n\npVecs = %.4f", mat.Formatted(pVecs, mat.Prefix("        ")))
-	fmt.Printf("\n\nqVecs = %.4f", mat.Formatted(qVecs, mat.Prefix("        ")))
+	fmt.Printf("\n\npVecs = %.4f", mat.Formatted(&pVecs, mat.Prefix("        ")))
+	fmt.Printf("\n\nqVecs = %.4f", mat.Formatted(&qVecs, mat.Prefix("        ")))
 
 	// Canonical Correlation Transforms. These can be useful as they represent
 	// the canonical variables as linear combinations of the original variables.
-	fmt.Printf("\n\nphiVs = %.4f", mat.Formatted(phiVs, mat.Prefix("        ")))
-	fmt.Printf("\n\npsiVs = %.4f", mat.Formatted(psiVs, mat.Prefix("        ")))
+	fmt.Printf("\n\nphiVs = %.4f", mat.Formatted(&phiVs, mat.Prefix("        ")))
+	fmt.Printf("\n\npsiVs = %.4f", mat.Formatted(&psiVs, mat.Prefix("        ")))
 
 	// Output:
 	// corRaw = ⎡-0.2192   0.3527   0.5828  -0.3883⎤
diff --git a/stat/cca_test.go b/stat/cca_test.go
index d3e6d62b..2429eb36 100644
--- a/stat/cca_test.go
+++ b/stat/cca_test.go
@@ -151,8 +151,8 @@ tests:
 	} {
 		var cc stat.CC
 		var corrs []float64
-		var pVecs, qVecs *mat.Dense
-		var phiVs, psiVs *mat.Dense
+		var pVecs, qVecs mat.Dense
+		var phiVs, psiVs mat.Dense
 		for j := 0; j < 2; j++ {
 			err := cc.CanonicalCorrelations(test.xdata, test.ydata, test.weights)
 			if err != nil {
@@ -161,30 +161,30 @@ tests:
 			}
 
 			corrs = cc.CorrsTo(corrs)
-			pVecs = cc.LeftTo(pVecs, true)
-			qVecs = cc.RightTo(qVecs, true)
-			phiVs = cc.LeftTo(phiVs, false)
-			psiVs = cc.RightTo(psiVs, false)
+			cc.LeftTo(&pVecs, true)
+			cc.RightTo(&qVecs, true)
+			cc.LeftTo(&phiVs, false)
+			cc.RightTo(&psiVs, false)
 
 			if !floats.EqualApprox(corrs, test.wantCorrs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected variance result got:%v, want:%v",
 					i, j, corrs, test.wantCorrs)
 			}
-			if !mat.EqualApprox(pVecs, test.wantpVecs, test.epsilon) {
+			if !mat.EqualApprox(&pVecs, test.wantpVecs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(pVecs), mat.Formatted(test.wantpVecs))
+					i, j, mat.Formatted(&pVecs), mat.Formatted(test.wantpVecs))
 			}
-			if !mat.EqualApprox(qVecs, test.wantqVecs, test.epsilon) {
+			if !mat.EqualApprox(&qVecs, test.wantqVecs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(qVecs), mat.Formatted(test.wantqVecs))
+					i, j, mat.Formatted(&qVecs), mat.Formatted(test.wantqVecs))
 			}
-			if !mat.EqualApprox(phiVs, test.wantphiVs, test.epsilon) {
+			if !mat.EqualApprox(&phiVs, test.wantphiVs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(phiVs), mat.Formatted(test.wantphiVs))
+					i, j, mat.Formatted(&phiVs), mat.Formatted(test.wantphiVs))
 			}
-			if !mat.EqualApprox(psiVs, test.wantpsiVs, test.epsilon) {
+			if !mat.EqualApprox(&psiVs, test.wantpsiVs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(psiVs), mat.Formatted(test.wantpsiVs))
+					i, j, mat.Formatted(&psiVs), mat.Formatted(test.wantpsiVs))
 			}
 		}
 	}
diff --git a/stat/combin/combin_test.go b/stat/combin/combin_test.go
index 50869e97..807e6f8e 100644
--- a/stat/combin/combin_test.go
+++ b/stat/combin/combin_test.go
@@ -298,6 +298,38 @@ func TestCartesian(t *testing.T) {
 	}
 }
 
+func TestNumCartesianProducts(t *testing.T) {
+	want := 6
+	got := Card([]int{1, 2, 3})
+	if want != got {
+		t.Errorf("number of cartesian products mismatch.\nwant:\n%v\ngot:\n%v", want, got)
+	}
+}
+
+func TestCartesianGenerator(t *testing.T) {
+	want := [][]int{
+		{0, 0, 0},
+		{0, 0, 1},
+		{0, 0, 2},
+		{0, 1, 0},
+		{0, 1, 1},
+		{0, 1, 2},
+	}
+	gen := NewCartesianGenerator([]int{1, 2, 3})
+	iterations := 0
+	for gen.Next() {
+		got := gen.Product()
+		if !reflect.DeepEqual(got, want[iterations]) {
+			t.Errorf("Cartesian product does not match. want: %v got: %v", want[iterations], got)
+		}
+		iterations++
+	}
+
+	if iterations != len(want) {
+		t.Errorf("Number of products does not match. want: %v got: %v", len(want), iterations)
+	}
+}
+
 func TestPermutationIndex(t *testing.T) {
 	for cas, s := range []struct {
 		n, k int
diff --git a/stat/combin/combinations_example_test.go b/stat/combin/combinations_example_test.go
index ff0cfeaa..4cba2ace 100644
--- a/stat/combin/combinations_example_test.go
+++ b/stat/combin/combinations_example_test.go
@@ -10,6 +10,49 @@ import (
 	"gonum.org/v1/gonum/stat/combin"
 )
 
+func ExampleCartesian() {
+	fmt.Println("Generate cartesian products for given lengths:")
+	lens := []int{1, 2, 3}
+	list := combin.Cartesian(lens)
+	for i, v := range list {
+		fmt.Println(i, v)
+	}
+	// This is easy, but the number of combinations  can be very large,
+	// and generating all at once can use a lot of memory.
+	// For big data sets, consider using CartesianGenerator instead.
+
+	// Output:
+	// Generate cartesian products for given lengths:
+	// 0 [0 0 0]
+	// 1 [0 0 1]
+	// 2 [0 0 2]
+	// 3 [0 1 0]
+	// 4 [0 1 1]
+	// 5 [0 1 2]
+}
+
+func ExampleCartesianGenerator() {
+	fmt.Println("Generate products for given lengths:")
+	lens := []int{1, 2, 3}
+	gen := combin.NewCartesianGenerator(lens)
+
+	// Now loop over all products.
+	var i int
+	for gen.Next() {
+		fmt.Println(i, gen.Product())
+		i++
+	}
+
+	// Output:
+	// Generate products for given lengths:
+	// 0 [0 0 0]
+	// 1 [0 0 1]
+	// 2 [0 0 2]
+	// 3 [0 1 0]
+	// 4 [0 1 1]
+	// 5 [0 1 2]
+}
+
 func ExampleCombinations() {
 	// combin provides several ways to work with the combinations of
 	// different objects. Combinations generates them directly.
diff --git a/stat/mds/mds_test.go b/stat/mds/mds_test.go
index 46da9846..8a9c8fd3 100644
--- a/stat/mds/mds_test.go
+++ b/stat/mds/mds_test.go
@@ -112,7 +112,8 @@ var torgersonScalingTests = []struct {
 func TestTorgersonScaling(t *testing.T) {
 	for i, test := range torgersonScalingTests {
 		_, c := test.dis.Dims()
-		gotK, got, gotVals := TorgersonScaling(nil, make([]float64, c), test.dis)
+		var got mat.Dense
+		gotK, gotVals := TorgersonScaling(&got, make([]float64, c), test.dis)
 		if gotK == 0 {
 			t.Error("unexpected scaling failure")
 			continue
@@ -120,9 +121,9 @@ func TestTorgersonScaling(t *testing.T) {
 		if gotK != test.wantK {
 			t.Errorf("unexpected k for test %d: got:%d want:%d", i, gotK, test.wantK)
 		}
-		if !mat.EqualApprox(colAbs{got}, colAbs{test.want}, 1e-5) {
+		if !mat.EqualApprox(colAbs{&got}, colAbs{test.want}, 1e-5) {
 			t.Errorf("unexpected result for test %d:\ngot:\n%.5f\nwant:\n%.5f",
-				i, mat.Formatted(got), mat.Formatted(test.want))
+				i, mat.Formatted(&got), mat.Formatted(test.want))
 		}
 		if !floats.EqualApprox(gotVals, test.wantVals, 1e-10) {
 			t.Errorf("unexpected Eigenvalues for test %d:\ngot: %.12e\nwant:%.12e", i, gotVals, test.wantVals)
diff --git a/stat/pca_example_test.go b/stat/pca_example_test.go
index 0ea7b50d..05ff56e2 100644
--- a/stat/pca_example_test.go
+++ b/stat/pca_example_test.go
@@ -40,7 +40,9 @@ func ExamplePC() {
 	// Project the data onto the first 2 principal components.
 	k := 2
 	var proj mat.Dense
-	proj.Mul(iris, pc.VectorsTo(nil).Slice(0, d, 0, k))
+	var vec mat.Dense
+	pc.VectorsTo(&vec)
+	proj.Mul(iris, vec.Slice(0, d, 0, k))
 
 	fmt.Printf("proj = %.4f", mat.Formatted(&proj, mat.Prefix("       ")))
 
diff --git a/stat/pca_test.go b/stat/pca_test.go
index 18356270..6eda9d50 100644
--- a/stat/pca_test.go
+++ b/stat/pca_test.go
@@ -151,11 +151,11 @@ tests:
 		},
 	} {
 		var pc PC
-		var vecs *mat.Dense
+		vecs := &mat.Dense{}
 		var vars []float64
 		for j := 0; j < 2; j++ {
 			ok := pc.PrincipalComponents(test.data, test.weights)
-			vecs = pc.VectorsTo(vecs)
+			pc.VectorsTo(vecs)
 			vars = pc.VarsTo(vars)
 			if !ok {
 				t.Errorf("unexpected SVD failure for test %d use %d", i, j)