native: add Dormhr

native/dormhr.go

@@ -0,0 +1,117 @@
+// 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 "github.com/gonum/blas"
+
+// Dormhr multiplies an m×n general matrix C with an nq×nq orthogonal matrix Q
+//  Q * C,    if side == blas.Left and trans == blas.NoTrans,
+//  Q^T * C,  if side == blas.Left and trans == blas.Trans,
+//  C * Q,    if side == blas.Right and trans == blas.NoTrans,
+//  C * Q^T,  if side == blas.Right and trans == blas.Trans,
+// where nq == m if side == blas.Left and nq == n if side == blas.Right.
+//
+// Q is defined implicitly as the product of ihi-ilo elementary reflectors, as
+// returned by Dgehrd:
+//  Q = H_{ilo} H_{ilo+1} ... H_{ihi-1}.
+// Q is equal to the identity matrix except in the submatrix
+// Q[ilo+1:ihi+1,ilo+1:ihi+1].
+//
+// ilo and ihi must have the same values as in the previous call of Dgehrd. It
+// must hold that
+//  0 <= ilo <= ihi < m,   if m > 0 and side == blas.Left,
+//  ilo = 0 and ihi = -1,  if m = 0 and side == blas.Left,
+//  0 <= ilo <= ihi < n,   if n > 0 and side == blas.Right,
+//  ilo = 0 and ihi = -1,  if n = 0 and side == blas.Right.
+//
+// a and lda represent an m×m matrix if side == blas.Left and an n×n matrix if
+// side == blas.Right. The matrix contains vectors which define the elementary
+// reflectors, as returned by Dgehrd.
+//
+// tau contains the scalar factors of the elementary reflectors, as returned by
+// Dgehrd. tau must have length m-1 if side == blas.Left and n-1 if side ==
+// blas.Right.
+//
+// c and ldc represent the m×n matrix C. On return, c is overwritten by the
+// product with Q.
+//
+// work must have length at least max(1,lwork), and lwork must be at least
+// max(1,n), if side == blas.Left, and max(1,m), if side == blas.Right. For
+// optimum performance lwork should be at least n*nb if side == blas.Left and
+// m*nb if side == blas.Right, where nb is the optimal block size. On return,
+// work[0] will contain the optimal value of lwork.
+//
+// If lwork == -1, instead of performing Dormhr, only the optimal value of lwork
+// will be stored in work[0].
+//
+// If any requirement on input sizes is not met, Dormhr will panic.
+//
+// Dormhr is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dormhr(side blas.Side, trans blas.Transpose, m, n, ilo, ihi int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
+	var (
+		nq int // The order of Q.
+		nw int // The minimum length of work.
+	)
+	switch side {
+	case blas.Left:
+		nq = m
+		nw = n
+	case blas.Right:
+		nq = n
+		nw = m
+	default:
+		panic(badSide)
+	}
+	checkMatrix(m, n, c, ldc)
+	checkMatrix(nq, nq, a, lda)
+	switch {
+	case trans != blas.NoTrans && trans != blas.Trans:
+		panic(badTrans)
+	case ilo < 0 || max(1, nq) <= ilo:
+		panic(badIlo)
+	case ihi < min(ilo, nq-1) || nq <= ihi:
+		panic(badIhi)
+	case nq > 0 && len(tau) != nq-1:
+		panic(badTau)
+	case lwork < max(1, nw) && lwork != -1:
+		panic(badWork)
+	case len(work) < max(1, lwork):
+		panic(shortWork)
+	}
+
+	nh := ihi - ilo
+	var nb int
+	if side == blas.Left {
+		opts := "LN"
+		if trans == blas.Trans {
+			opts = "LT"
+		}
+		nb = impl.Ilaenv(1, "DORMQR", opts, nh, n, nh, -1)
+	} else {
+		opts := "RN"
+		if trans == blas.Trans {
+			opts = "RT"
+		}
+		nb = impl.Ilaenv(1, "DORMQR", opts, m, nh, nh, -1)
+	}
+	lwkopt := max(1, nw) * nb
+	if lwork == -1 {
+		work[0] = float64(lwkopt)
+		return
+	}
+
+	if m == 0 || n == 0 || nh == 0 {
+		work[0] = 1
+		return
+	}
+	if side == blas.Left {
+		impl.Dormqr(side, trans, nh, n, nh, a[(ilo+1)*lda+ilo:], lda,
+			tau[ilo:ihi], c[(ilo+1)*ldc:], ldc, work, lwork)
+	} else {
+		impl.Dormqr(side, trans, m, nh, nh, a[(ilo+1)*lda+ilo:], lda,
+			tau[ilo:ihi], c[ilo+1:], ldc, work, lwork)
+	}
+	work[0] = float64(lwkopt)
+}