native: add Dlaqr1

native/dlaqr1.go

@@ -0,0 +1,57 @@
+// Copyright ©2016 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 native
+
+import "math"
+
+// Dlaqr1 sets v to a scalar multiple of the first column of the product
+//  (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
+// where H is a 2×2 or 3×3 matrix, I is the identity matrix of the same size,
+// and i is the imaginary unit. Scaling is done to avoid overflows and most
+// underflows.
+//
+// n is the order of H and must be either 2 or 3. It must hold that either sr1 =
+// sr2 and si1 = -si2, or si1 = si2 = 0. The length of v must be equal to n. If
+// any of these conditions is not met, Dlaqr1 will panic.
+//
+// Dlaqr1 is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64) {
+	if n != 2 && n != 3 {
+		panic(badDims)
+	}
+	checkMatrix(n, n, h, ldh)
+	if len(v) != n {
+		panic(badSlice)
+	}
+	if !((sr1 == sr2 && si1 == -si2) || (si1 == 0 && si2 == 0)) {
+		panic(badShifts)
+	}
+
+	if n == 2 {
+		s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh])
+		if s == 0 {
+			v[0] = 0
+			v[1] = 0
+		} else {
+			h21s := h[ldh] / s
+			v[0] = h21s*h[1] + (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s)
+			v[1] = h21s * (h[0] + h[ldh+1] - sr1 - sr2)
+		}
+		return
+	}
+
+	s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh]) + math.Abs(h[2*ldh])
+	if s == 0 {
+		v[0] = 0
+		v[1] = 0
+		v[2] = 0
+	} else {
+		h21s := h[ldh] / s
+		h31s := h[2*ldh] / s
+		v[0] = (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s) + h[1]*h21s + h[2]*h31s
+		v[1] = h21s*(h[0]+h[ldh+1]-sr1-sr2) + h[ldh+2]*h31s
+		v[2] = h31s*(h[0]+h[2*ldh+2]-sr1-sr2) + h21s*h[2*ldh+1]
+	}
+}

native/general.go

@@ -32,6 +32,7 @@ const (
 	badNorm         = "lapack: bad norm"
 	badPivot        = "lapack: bad pivot"
 	badS            = "lapack: s has insufficient length"
+	badShifts       = "lapack: bad shifts"
 	badSide         = "lapack: bad side"
 	badSlice        = "lapack: bad input slice length"
 	badStore        = "lapack: bad store"