Repository: gonum/gonum
Pull request #1115 (2019-10-07 23:39:54)

graph/path: add nodes to shortest path tree lazily
This is an alternative to #1113 that gives the full benefits of the change to lazy addition to the shortest path tree during `AStar` execution, while also retaining the implied semantics of `graph.Graph` (that the graph does not change under the callee's feet).

This reduces total volume of allocations, particularly in the case where there is a heuristic provided at the cost of allocating more often in general. In the same situation, the time taken is dramatically improved. It also reduces the requirement of the input graph's method set (to match now what `Dijkstra` needs).

Please take a look.

Closes #1113.

/cc @phyrwork

<details>
<table class='benchstat oldnew'>
<tr class='configs'><th><th>old<th>new


<tbody>
<tr><th><th colspan='2' class='metric'>time/op<th>delta
<tr class='better'><td>AStarUndirected/GNP_Undirected_10_tenth-8<td>3.86µs ± 3%<td>3.44µs ± 1%<td class='delta'>−10.81%<td class='note'>(p=0.000 n=9&#43;8)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_tenth-8<td>215µs ± 0%<td>217µs ± 0%<td class='delta'>&#43;1.14%<td class='note'>(p=0.000 n=9&#43;7)
<tr class='unchanged'><td>AStarUndirected/GNP_Undirected_1000_tenth-8<td>10.9ms ± 8%<td>10.8ms ± 4%<td class='nodelta'>~<td class='note'>(p=0.696 n=10&#43;8)
<tr class='better'><td>AStarUndirected/GNP_Undirected_10_half-8<td>9.60µs ± 0%<td>9.57µs ± 1%<td class='delta'>−0.35%<td class='note'>(p=0.041 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_half-8<td>669µs ± 1%<td>676µs ± 1%<td class='delta'>&#43;1.01%<td class='note'>(p=0.001 n=10&#43;10)
<tr class='unchanged'><td>AStarUndirected/GNP_Undirected_1000_half-8<td>28.6ms ± 8%<td>28.4ms ± 6%<td class='nodelta'>~<td class='note'>(p=0.529 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_2_2-8<td>52.5µs ± 1%<td>60.6µs ± 1%<td class='delta'>&#43;15.37%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_10_2_2_2_heuristic-8<td>15.9µs ± 0%<td>11.2µs ± 0%<td class='delta'>−29.79%<td class='note'>(p=0.000 n=9&#43;8)
<tr class='better'><td>AStarUndirected/NSW_Undirected_10_2_5_2-8<td>77.9µs ± 1%<td>76.7µs ± 1%<td class='delta'>−1.61%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_10_2_5_2_heuristic-8<td>16.3µs ± 1%<td>11.7µs ± 2%<td class='delta'>−28.30%<td class='note'>(p=0.000 n=9&#43;9)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_10_2-8<td>1.78ms ± 3%<td>1.34ms ± 3%<td class='delta'>−24.91%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_10_2_heuristic-8<td>718µs ± 3%<td>53µs ± 1%<td class='delta'>−92.65%<td class='note'>(p=0.000 n=8&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_20_2-8<td>4.37ms ± 5%<td>3.99ms ± 5%<td class='delta'>−8.59%<td class='note'>(p=0.000 n=9&#43;9)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_20_2_heuristic-8<td>759µs ± 1%<td>121µs ±33%<td class='delta'>−84.08%<td class='note'>(p=0.000 n=10&#43;10)
<tr><td>&nbsp;
</tbody>

<tbody>
<tr><th><th colspan='2' class='metric'>alloc/op<th>delta
<tr class='better'><td>AStarUndirected/GNP_Undirected_10_tenth-8<td>1.65kB ± 0%<td>1.60kB ± 0%<td class='delta'>−3.15%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_tenth-8<td>51.6kB ± 0%<td>58.7kB ± 0%<td class='delta'>&#43;13.72%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_1000_tenth-8<td>1.82MB ± 8%<td>1.89MB ± 5%<td class='delta'>&#43;4.22%<td class='note'>(p=0.043 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_10_half-8<td>3.14kB ± 0%<td>3.89kB ± 0%<td class='delta'>&#43;23.91%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_half-8<td>127kB ± 1%<td>136kB ± 1%<td class='delta'>&#43;6.82%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='unchanged'><td>AStarUndirected/GNP_Undirected_1000_half-8<td>4.17MB ± 8%<td>4.25MB ± 6%<td class='nodelta'>~<td class='note'>(p=0.393 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_2_2-8<td>20.7kB ± 0%<td>26.3kB ± 0%<td class='delta'>&#43;27.10%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_10_2_2_2_heuristic-8<td>10.8kB ± 0%<td>7.3kB ± 0%<td class='delta'>−32.14%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_5_2-8<td>28.2kB ± 0%<td>33.6kB ± 0%<td class='delta'>&#43;18.97%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_10_2_5_2_heuristic-8<td>10.9kB ± 0%<td>7.4kB ± 0%<td class='delta'>−31.84%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_10_2-8<td>1.11MB ± 2%<td>0.81MB ± 3%<td class='delta'>−26.77%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_10_2_heuristic-8<td>668kB ± 0%<td>32kB ± 0%<td class='delta'>−95.21%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_100_5_20_2-8<td>2.00MB ± 4%<td>2.12MB ± 3%<td class='delta'>&#43;6.05%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='better'><td>AStarUndirected/NSW_Undirected_100_5_20_2_heuristic-8<td>687kB ± 0%<td>64kB ± 0%<td class='delta'>−90.64%<td class='note'>(p=0.000 n=9&#43;10)
<tr><td>&nbsp;
</tbody>

<tbody>
<tr><th><th colspan='2' class='metric'>allocs/op<th>delta
<tr class='worse'><td>AStarUndirected/GNP_Undirected_10_tenth-8<td>33.0 ± 0%<td>38.0 ± 0%<td class='delta'>&#43;15.15%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_tenth-8<td>621 ± 0%<td>646 ± 0%<td class='delta'>&#43;4.08%<td class='note'>(p=0.000 n=10&#43;8)
<tr class='unchanged'><td>AStarUndirected/GNP_Undirected_1000_tenth-8<td>20.8k ± 7%<td>20.7k ± 5%<td class='nodelta'>~<td class='note'>(p=0.684 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_10_half-8<td>57.0 ± 0%<td>66.0 ± 0%<td class='delta'>&#43;15.79%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/GNP_Undirected_100_half-8<td>1.47k ± 1%<td>1.51k ± 1%<td class='delta'>&#43;2.81%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='unchanged'><td>AStarUndirected/GNP_Undirected_1000_half-8<td>55.9k ± 8%<td>56.0k ± 6%<td class='nodelta'>~<td class='note'>(p=0.971 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_2_2-8<td>149 ± 0%<td>188 ± 1%<td class='delta'>&#43;25.91%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_2_2_heuristic-8<td>47.0 ± 0%<td>56.0 ± 0%<td class='delta'>&#43;19.15%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_5_2-8<td>212 ± 0%<td>230 ± 0%<td class='delta'>&#43;8.53%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_10_2_5_2_heuristic-8<td>48.0 ± 0%<td>58.0 ± 0%<td class='delta'>&#43;20.83%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_100_5_10_2-8<td>2.80k ± 4%<td>3.06k ± 2%<td class='delta'>&#43;9.55%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_100_5_10_2_heuristic-8<td>141 ± 0%<td>163 ± 0%<td class='delta'>&#43;15.60%<td class='note'>(p=0.000 n=10&#43;10)
<tr class='unchanged'><td>AStarUndirected/NSW_Undirected_100_5_20_2-8<td>7.99k ± 6%<td>7.75k ± 3%<td class='nodelta'>~<td class='note'>(p=0.075 n=10&#43;10)
<tr class='worse'><td>AStarUndirected/NSW_Undirected_100_5_20_2_heuristic-8<td>242 ± 0%<td>275 ± 0%<td class='delta'>&#43;13.64%<td class='note'>(p=0.000 n=8&#43;10)
<tr><td>&nbsp;
</tbody>

</table>
</details>

<!--
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 _____."
-->



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 194f39bf..b407b2ae 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -65,6 +65,7 @@ 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>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 1b16a4ff..1b6e3414 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -71,6 +71,7 @@ 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>
diff --git a/README.md b/README.md
index 81c25b2f..9b49823d 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@ 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 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).
 
 
 ## Build tags
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/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/simpsons.go b/integrate/simpsons.go
new file mode 100644
index 00000000..773845bd
--- /dev/null
+++ b/integrate/simpsons.go
@@ -0,0 +1,90 @@
+// 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"
+
+// Simpson's method approximates the integral of a function f(x)
+// by means of 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_{2i} + {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 * f_{N-2} + {a_1}_i * f_{N-1} + {a_2}_i * f_{N}
+// where the coefficients are:
+//  {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 function of f if N is an even number of subintervals.
+// If the number of subintervals are odd, the subintervals up i=0..n-2 are given by the above
+// and the approximations over the second to last and last subintervals are 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})
+// More information is available at:
+// https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
+//
+// The (x,f) input data points must be sorted along x. One can use stat.SortWeighted to do that.
+// The x and f slices must be of equal length and the length must be greater than two.
+func Simpsons(x, f []float64) float64 {
+	switch {
+	case len(x) != len(f):
+		panic("integrate: slice length mismatch")
+	case len(x) <= 2:
+		panic("integrate: input data too small")
+	case !sort.Float64sAreSorted(x):
+		panic("integrate: must be sorted")
+	}
+
+	var integral float64
+	for i := 1; i < len(x)-1; i += 2 {
+		if x[i] == x[i-1] || x[i] == x[i+1] {
+			panic("integrate: repeated abscissa are not allowed")
+		}
+		h0 := x[i] - x[i-1]
+		h0p2 := h0 * h0
+		h0p3 := h0 * h0 * h0
+		h1 := x[i+1] - x[i]
+		h1p2 := h1 * h1
+		h1p3 := h1 * h1 * h1
+		hph := h0 + h1
+		f0 := f[i-1]
+		f1 := f[i]
+		f2 := f[i+1]
+		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 * f0
+		integral += a1 * f1
+		integral += a2 * f2
+	}
+
+	if len(x)%2 == 0 {
+		h0 := x[len(x)-2] - x[len(x)-3]
+		h1 := x[len(x)-1] - x[len(x)-2]
+		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[len(x)-3]
+		integral += a1 * f[len(x)-2]
+		integral += a2 * f[len(x)-1]
+	}
+
+	return integral
+}
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/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/simpsons_test.go b/integrate/simpsons_test.go
new file mode 100644
index 00000000..cdaff0c2
--- /dev/null
+++ b/integrate/simpsons_test.go
@@ -0,0 +1,172 @@
+// 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"
+)
+
+func TestSimpsons(t *testing.T) {
+	x := floats.Span(make([]float64, 1e6), 0, 1)
+	for i, test := range []struct {
+		x    []float64
+		f    func(x float64) float64
+		want float64
+		tol  float64
+	}{
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 1),
+			f:    func(x float64) float64 { return math.Pi },
+			want: math.Pi,
+			tol:  1e-10,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 1),
+			f:    func(x float64) float64 { return 1.0 },
+			want: 1,
+			tol:  1e-10,
+		},
+		{
+			x:    floats.Span(make([]float64, 3), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 10), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return x },
+			want: 0.5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), -1, 1),
+			f:    func(x float64) float64 { return x },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return x + 10 },
+			want: 10.5,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return 3*x*x + 10 },
+			want: 11,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return math.Exp(x) },
+			want: 1.7182818284591876,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, math.Pi),
+			f:    func(x float64) float64 { return math.Cos(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
+			f:    func(x float64) float64 { return math.Cos(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6*10), 0, math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6*10), 0, 0.5*math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 1,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x: join(floats.Span(make([]float64, 3), 0, math.Pi/3),
+				[]float64{4 * math.Pi / 10},
+				floats.Span(make([]float64, 4), 3*math.Pi/4, math.Pi)),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-2,
+		},
+		{
+			x: join(floats.Span(make([]float64, 30), 0, 5*math.Pi/16),
+				floats.Span(make([]float64, 100), 3*math.Pi/8, math.Pi)),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-4,
+		},
+		{
+			x: join(floats.Span(make([]float64, 1e5), 0, 15),
+				floats.Span(make([]float64, 1e5), 23, 40),
+				floats.Span(make([]float64, 1e5), 50, 80),
+				floats.Span(make([]float64, 1e5), 90, 100)),
+			f:    func(x float64) float64 { return 2 * x },
+			want: 10000,
+			tol:  1e-9,
+		},
+		{
+			x:    []float64{0, 1, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+		{
+			x:    []float64{0, 0.2, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+		{
+			x:    []float64{0, 1.89, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+	} {
+		y := make([]float64, len(test.x))
+		for i, v := range test.x {
+			y[i] = test.f(v)
+		}
+		v := Simpsons(test.x, y)
+		if !floats.EqualWithinAbs(v, test.want, test.tol) {
+			t.Errorf("test #%d: got=%v want=%f\n", i, v, float64(test.want))
+		}
+	}
+}
+
+func join(slices ...[]float64) []float64 {
+	var c []float64
+	for _, s := range slices {
+		c = append(c, s...)
+	}
+	return c
+}
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/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)