Golang blas64.Implementation函数代码示例

// Drscl multiplies the vector x by 1/a being careful to avoid overflow or
// underflow where possible.
func (impl Implementation) Drscl(n int, a float64, x []float64, incX int) {
	checkVector(n, x, incX)
	bi := blas64.Implementation()
	cden := a
	cnum := 1.0
	smlnum := dlamchS
	bignum := 1 / smlnum
	for {
		cden1 := cden * smlnum
		cnum1 := cnum / bignum
		var mul float64
		var done bool
		switch {
		case cnum != 0 && math.Abs(cden1) > math.Abs(cnum):
			mul = smlnum
			done = false
			cden = cden1
		case math.Abs(cnum1) > math.Abs(cden):
			mul = bignum
			done = false
			cnum = cnum1
		default:
			mul = cnum / cden
			done = true
		}
		bi.Dscal(n, mul, x, incX)
		if done {
			break
		}
	}
}

// Dgetrs solves a system of equations using an LU factorization.
// The system of equations solved is
//  A * X = B if trans == blas.Trans
//  A^T * X = B if trans == blas.NoTrans
// A is a general n×n matrix with stride lda. B is a general matrix of size n×nrhs.
//
// On entry b contains the elements of the matrix B. On exit, b contains the
// elements of X, the solution to the system of equations.
//
// a and ipiv contain the LU factorization of A and the permutation indices as
// computed by Dgetrf. ipiv is zero-indexed.
func (impl Implementation) Dgetrs(trans blas.Transpose, n, nrhs int, a []float64, lda int, ipiv []int, b []float64, ldb int) {
	checkMatrix(n, n, a, lda)
	checkMatrix(n, nrhs, b, ldb)
	if len(ipiv) < n {
		panic(badIpiv)
	}
	if n == 0 || nrhs == 0 {
		return
	}
	if trans != blas.Trans && trans != blas.NoTrans {
		panic(badTrans)
	}
	bi := blas64.Implementation()
	if trans == blas.NoTrans {
		// Solve A * X = B.
		impl.Dlaswp(nrhs, b, ldb, 0, n-1, ipiv, 1)
		// Solve L * X = B, updating b.
		bi.Dtrsm(blas.Left, blas.Lower, blas.NoTrans, blas.Unit,
			n, nrhs, 1, a, lda, b, ldb)
		// Solve U * X = B, updating b.
		bi.Dtrsm(blas.Left, blas.Upper, blas.NoTrans, blas.NonUnit,
			n, nrhs, 1, a, lda, b, ldb)
		return
	}
	// Solve A^T * X = B.
	// Solve U^T * X = B, updating b.
	bi.Dtrsm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit,
		n, nrhs, 1, a, lda, b, ldb)
	// Solve L^T * X = B, updating b.
	bi.Dtrsm(blas.Left, blas.Lower, blas.Trans, blas.Unit,
		n, nrhs, 1, a, lda, b, ldb)
	impl.Dlaswp(nrhs, b, ldb, 0, n-1, ipiv, -1)
}

// Dlarf applies an elementary reflector to a general rectangular matrix c.
// This computes
//  c = h * c if side == Left
//  c = c * h if side == right
// where
//  h = 1 - tau * v * v^T
// and c is an m * n matrix.
//
// work is temporary storage of length at least m if side == Left and at least
// n if side == Right. This function will panic if this length requirement is not met.
//
// Dlarf is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
	applyleft := side == blas.Left
	if (applyleft && len(work) < n) || (!applyleft && len(work) < m) {
		panic(badWork)
	}
	checkMatrix(m, n, c, ldc)

	// v has length m if applyleft and n otherwise.
	lenV := n
	if applyleft {
		lenV = m
	}

	checkVector(lenV, v, incv)

	lastv := 0 // last non-zero element of v
	lastc := 0 // last non-zero row/column of c
	if tau != 0 {
		var i int
		if applyleft {
			lastv = m - 1
		} else {
			lastv = n - 1
		}
		if incv > 0 {
			i = lastv * incv
		}

		// Look for the last non-zero row in v.
		for lastv >= 0 && v[i] == 0 {
			lastv--
			i -= incv
		}
		if applyleft {
			// Scan for the last non-zero column in C[0:lastv, :]
			lastc = impl.Iladlc(lastv+1, n, c, ldc)
		} else {
			// Scan for the last non-zero row in C[:, 0:lastv]
			lastc = impl.Iladlr(m, lastv+1, c, ldc)
		}
	}
	if lastv == -1 || lastc == -1 {
		return
	}
	// Sometimes 1-indexing is nicer ...
	bi := blas64.Implementation()
	if applyleft {
		// Form H * C
		// w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]^T * v[1:lastv+1,1]
		bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
		// c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]^T
		bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
		return
	}
	// Form C*H
	// w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
	bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
	// c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]^T
	bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
}

// Dpotrf computes the cholesky decomposition of the symmetric positive definite
// matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
// and a = U U^T is stored in place into a. If ul == blas.Lower, then a = L L^T
// is computed and stored in-place into a. If a is not positive definite, false
// is returned. This is the blocked version of the algorithm.
func (impl Implementation) Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
	bi := blas64.Implementation()
	if ul != blas.Upper && ul != blas.Lower {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if lda < n {
		panic(badLdA)
	}
	if n == 0 {
		return true
	}
	nb := impl.Ilaenv(1, "DPOTRF", string(ul), n, -1, -1, -1)
	if n <= nb {
		return impl.Dpotf2(ul, n, a, lda)
	}
	if ul == blas.Upper {
		for j := 0; j < n; j += nb {
			jb := min(nb, n-j)
			bi.Dsyrk(blas.Upper, blas.Trans, jb, j,
				-1, a[j:], lda,
				1, a[j*lda+j:], lda)
			ok = impl.Dpotf2(blas.Upper, jb, a[j*lda+j:], lda)
			if !ok {
				return ok
			}
			if j+jb < n {
				bi.Dgemm(blas.Trans, blas.NoTrans, jb, n-j-jb, j,
					-1, a[j:], lda, a[j+jb:], lda,
					1, a[j*lda+j+jb:], lda)
				bi.Dtrsm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit, jb, n-j-jb,
					1, a[j*lda+j:], lda,
					a[j*lda+j+jb:], lda)
			}
		}
		return true
	}
	for j := 0; j < n; j += nb {
		jb := min(nb, n-j)
		bi.Dsyrk(blas.Lower, blas.NoTrans, jb, j,
			-1, a[j*lda:], lda,
			1, a[j*lda+j:], lda)
		ok := impl.Dpotf2(blas.Lower, jb, a[j*lda+j:], lda)
		if !ok {
			return ok
		}
		if j+jb < n {
			bi.Dgemm(blas.NoTrans, blas.Trans, n-j-jb, jb, j,
				-1, a[(j+jb)*lda:], lda, a[j*lda:], lda,
				1, a[(j+jb)*lda+j:], lda)
			bi.Dtrsm(blas.Right, blas.Lower, blas.Trans, blas.NonUnit, n-j-jb, jb,
				1, a[j*lda+j:], lda,
				a[(j+jb)*lda+j:], lda)
		}
	}
	return true
}

// Dgecon estimates the reciprocal of the condition number of the n×n matrix A
// given the LU decomposition of the matrix. The condition number computed may
// be based on the 1-norm or the ∞-norm.
//
// The slice a contains the result of the LU decomposition of A as computed by Dgetrf.
//
// anorm is the corresponding 1-norm or ∞-norm of the original matrix A.
//
// work is a temporary data slice of length at least 4*n and Dgecon will panic otherwise.
//
// iwork is a temporary data slice of length at least n and Dgecon will panic otherwise.
func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
	checkMatrix(n, n, a, lda)
	if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum {
		panic(badNorm)
	}
	if len(work) < 4*n {
		panic(badWork)
	}
	if len(iwork) < n {
		panic(badWork)
	}

	if n == 0 {
		return 1
	} else if anorm == 0 {
		return 0
	}

	bi := blas64.Implementation()
	var rcond, ainvnm float64
	var kase int
	var normin bool
	isave := new([3]int)
	onenrm := norm == lapack.MaxColumnSum
	smlnum := dlamchS
	kase1 := 2
	if onenrm {
		kase1 = 1
	}
	for {
		ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
		if kase == 0 {
			if ainvnm != 0 {
				rcond = (1 / ainvnm) / anorm
			}
			return rcond
		}
		var sl, su float64
		if kase == kase1 {
			sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.Unit, normin, n, a, lda, work, work[2*n:])
			su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
		} else {
			su = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
			sl = impl.Dlatrs(blas.Lower, blas.Trans, blas.Unit, normin, n, a, lda, work, work[2*n:])
		}
		scale := sl * su
		normin = true
		if scale != 1 {
			ix := bi.Idamax(n, work, 1)
			if scale == 0 || scale < math.Abs(work[ix])*smlnum {
				return rcond
			}
			impl.Drscl(n, scale, work, 1)
		}
	}
}

// Dgebak updates an n×m matrix V as
//  V = P D V,        if side == blas.Right,
//  V = P D^{-1} V,   if side == blas.Left,
// where P and D are n×n permutation and scaling matrices, respectively,
// implicitly represented by job, scale, ilo and ihi as returned by Dgebal.
//
// Typically, columns of the matrix V contain the right or left (determined by
// side) eigenvectors of the balanced matrix output by Dgebal, and Dgebak forms
// the eigenvectors of the original matrix.
//
// Dgebak is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dgebak(job lapack.Job, side blas.Side, n, ilo, ihi int, scale []float64, m int, v []float64, ldv int) {
	switch job {
	default:
		panic(badJob)
	case lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale:
	}
	switch side {
	default:
		panic(badSide)
	case blas.Left, blas.Right:
	}
	checkMatrix(n, m, v, ldv)
	switch {
	case ilo < 0 || max(0, n-1) < ilo:
		panic(badIlo)
	case ihi < min(ilo, n-1) || n <= ihi:
		panic(badIhi)
	}

	// Quick return if possible.
	if n == 0 || m == 0 || job == lapack.None {
		return
	}

	bi := blas64.Implementation()
	if ilo != ihi && job != lapack.Permute {
		// Backward balance.
		if side == blas.Right {
			for i := ilo; i <= ihi; i++ {
				bi.Dscal(m, scale[i], v[i*ldv:], 1)
			}
		} else {
			for i := ilo; i <= ihi; i++ {
				bi.Dscal(m, 1/scale[i], v[i*ldv:], 1)
			}
		}
	}
	if job == lapack.Scale {
		return
	}
	// Backward permutation.
	for i := ilo - 1; i >= 0; i-- {
		k := int(scale[i])
		if k == i {
			continue
		}
		bi.Dswap(m, v[i*ldv:], 1, v[k*ldv:], 1)
	}
	for i := ihi + 1; i < n; i++ {
		k := int(scale[i])
		if k == i {
			continue
		}
		bi.Dswap(m, v[i*ldv:], 1, v[k*ldv:], 1)
	}
}

// Dpotf2 computes the cholesky decomposition of the symmetric positive definite
// matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
// and a = U^T U is stored in place into a. If ul == blas.Lower, then a = L L^T
// is computed and stored in-place into a. If a is not positive definite, false
// is returned. This is the unblocked version of the algorithm.
func (Implementation) Dpotf2(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
	if ul != blas.Upper && ul != blas.Lower {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if lda < n {
		panic(badLdA)
	}
	if n == 0 {
		return true
	}
	bi := blas64.Implementation()
	if ul == blas.Upper {
		for j := 0; j < n; j++ {
			ajj := a[j*lda+j]
			if j != 0 {
				ajj -= bi.Ddot(j, a[j:], lda, a[j:], lda)
			}
			if ajj <= 0 || math.IsNaN(ajj) {
				a[j*lda+j] = ajj
				return false
			}
			ajj = math.Sqrt(ajj)
			a[j*lda+j] = ajj
			if j < n-1 {
				bi.Dgemv(blas.Trans, j, n-j-1,
					-1, a[j+1:], lda, a[j:], lda,
					1, a[j*lda+j+1:], 1)
				bi.Dscal(n-j-1, 1/ajj, a[j*lda+j+1:], 1)
			}
		}
		return true
	}
	for j := 0; j < n; j++ {
		ajj := a[j*lda+j]
		if j != 0 {
			ajj -= bi.Ddot(j, a[j*lda:], 1, a[j*lda:], 1)
		}
		if ajj <= 0 || math.IsNaN(ajj) {
			a[j*lda+j] = ajj
			return false
		}
		ajj = math.Sqrt(ajj)
		a[j*lda+j] = ajj
		if j < n-1 {
			bi.Dgemv(blas.NoTrans, n-j-1, j,
				-1, a[(j+1)*lda:], lda, a[j*lda:], 1,
				1, a[(j+1)*lda+j:], lda)
			bi.Dscal(n-j-1, 1/ajj, a[(j+1)*lda+j:], lda)
		}
	}
	return true
}

// Dpocon estimates the reciprocal of the condition number of a positive-definite
// matrix A given the Cholesky decmposition of A. The condition number computed
// is based on the 1-norm and the ∞-norm.
//
// anorm is the 1-norm and the ∞-norm of the original matrix A.
//
// work is a temporary data slice of length at least 3*n and Dpocon will panic otherwise.
//
// iwork is a temporary data slice of length at least n and Dpocon will panic otherwise.
func (impl Implementation) Dpocon(uplo blas.Uplo, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
	checkMatrix(n, n, a, lda)
	if uplo != blas.Upper && uplo != blas.Lower {
		panic(badUplo)
	}
	if len(work) < 3*n {
		panic(badWork)
	}
	if len(iwork) < n {
		panic(badWork)
	}
	var rcond float64
	if n == 0 {
		return 1
	}
	if anorm == 0 {
		return rcond
	}

	bi := blas64.Implementation()
	var ainvnm float64
	smlnum := dlamchS
	upper := uplo == blas.Upper
	var kase int
	var normin bool
	isave := new([3]int)
	var sl, su float64
	for {
		ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
		if kase == 0 {
			if ainvnm != 0 {
				rcond = (1 / ainvnm) / anorm
			}
			return rcond
		}
		if upper {
			sl = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
			normin = true
			su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
		} else {
			sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
			normin = true
			su = impl.Dlatrs(blas.Lower, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
		}
		scale := sl * su
		if scale != 1 {
			ix := bi.Idamax(n, work, 1)
			if scale == 0 || scale < math.Abs(work[ix])*smlnum {
				return rcond
			}
			impl.Drscl(n, scale, work, 1)
		}
	}
}

func DpotrfTest(t *testing.T, impl Dpotrfer) {
	rnd := rand.New(rand.NewSource(1))
	bi := blas64.Implementation()
	for i, test := range []struct {
		n int
	}{
		{n: 10},
		{n: 30},
		{n: 63},
		{n: 65},
		{n: 128},
		{n: 1000},
	} {
		n := test.n
		// Construct a positive-definite symmetric matrix
		base := make([]float64, n*n)
		for i := range base {
			base[i] = rnd.Float64()
		}
		a := make([]float64, len(base))
		bi.Dgemm(blas.Trans, blas.NoTrans, n, n, n, 1, base, n, base, n, 0, a, n)

		aCopy := make([]float64, len(a))
		copy(aCopy, a)

		// Test with Upper
		impl.Dpotrf(blas.Upper, n, a, n)

		// zero all the other elements
		for i := 0; i < n; i++ {
			for j := 0; j < i; j++ {
				a[i*n+j] = 0
			}
		}
		// Multiply u^T * u
		ans := make([]float64, len(a))
		bi.Dsyrk(blas.Upper, blas.Trans, n, n, 1, a, n, 0, ans, n)

		match := true
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				if !floats.EqualWithinAbsOrRel(ans[i*n+j], aCopy[i*n+j], 1e-14, 1e-14) {
					match = false
				}
			}
		}
		if !match {
			//fmt.Println(aCopy)
			//fmt.Println(ans)
			t.Errorf("Case %v: Mismatch for upper", i)
		}
	}
}

// Dlaswp swaps the rows k1 to k2 of a according to the indices in ipiv.
// a is a matrix with n columns and stride lda. incX is the increment for ipiv.
// k1 and k2 are zero-indexed. If incX is negative, then loops from k2 to k1
//
// Dlaswp is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlaswp(n int, a []float64, lda, k1, k2 int, ipiv []int, incX int) {
	if incX != 1 && incX != -1 {
		panic(absIncNotOne)
	}
	bi := blas64.Implementation()
	if incX == 1 {
		for k := k1; k <= k2; k++ {
			bi.Dswap(n, a[k*lda:], 1, a[ipiv[k]*lda:], 1)
		}
		return
	}
	for k := k2; k >= k1; k-- {
		bi.Dswap(n, a[k*lda:], 1, a[ipiv[k]*lda:], 1)
	}
}

// Dtrtri computes the inverse of a triangular matrix, storing the result in place
// into a. This is the BLAS level 3 version of the algorithm which builds upon
// Dtrti2 to operate on matrix blocks instead of only individual columns.
//
// Dtrtri will not perform the inversion if the matrix is singular, and returns
// a boolean indicating whether the inversion was successful.
func (impl Implementation) Dtrtri(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) (ok bool) {
	checkMatrix(n, n, a, lda)
	if uplo != blas.Upper && uplo != blas.Lower {
		panic(badUplo)
	}
	if diag != blas.NonUnit && diag != blas.Unit {
		panic(badDiag)
	}
	if n == 0 {
		return false
	}
	nonUnit := diag == blas.NonUnit
	if nonUnit {
		for i := 0; i < n; i++ {
			if a[i*lda+i] == 0 {
				return false
			}
		}
	}

	bi := blas64.Implementation()

	nb := impl.Ilaenv(1, "DTRTRI", "UD", n, -1, -1, -1)
	if nb <= 1 || nb > n {
		impl.Dtrti2(uplo, diag, n, a, lda)
		return true
	}
	if uplo == blas.Upper {
		for j := 0; j < n; j += nb {
			jb := min(nb, n-j)
			bi.Dtrmm(blas.Left, blas.Upper, blas.NoTrans, diag, j, jb, 1, a, lda, a[j:], lda)
			bi.Dtrsm(blas.Right, blas.Upper, blas.NoTrans, diag, j, jb, -1, a[j*lda+j:], lda, a[j:], lda)
			impl.Dtrti2(blas.Upper, diag, jb, a[j*lda+j:], lda)
		}
		return true
	}
	nn := ((n - 1) / nb) * nb
	for j := nn; j >= 0; j -= nb {
		jb := min(nb, n-j)
		if j+jb <= n-1 {
			bi.Dtrmm(blas.Left, blas.Lower, blas.NoTrans, diag, n-j-jb, jb, 1, a[(j+jb)*lda+j+jb:], lda, a[(j+jb)*lda+j:], lda)
			bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, diag, n-j-jb, jb, -1, a[j*lda+j:], lda, a[(j+jb)*lda+j:], lda)
		}
		impl.Dtrti2(blas.Lower, diag, jb, a[j*lda+j:], lda)
	}
	return true
}

// Dorg2r generates an m×n matrix Q with orthonormal columns defined by the
// product of elementary reflectors as computed by Dgeqrf.
//  Q = H_0 * H_1 * ... * H_{k-1}
// len(tau) >= k, 0 <= k <= n, 0 <= n <= m, len(work) >= n.
// Dorg2r will panic if these conditions are not met.
//
// Dorg2r is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dorg2r(m, n, k int, a []float64, lda int, tau []float64, work []float64) {
	checkMatrix(m, n, a, lda)
	if len(tau) < k {
		panic(badTau)
	}
	if len(work) < n {
		panic(badWork)
	}
	if k > n {
		panic(kGTN)
	}
	if n > m {
		panic(mLTN)
	}
	if len(work) < n {
		panic(badWork)
	}
	if n == 0 {
		return
	}
	bi := blas64.Implementation()
	// Initialize columns k+1:n to columns of the unit matrix.
	for l := 0; l < m; l++ {
		for j := k; j < n; j++ {
			a[l*lda+j] = 0
		}
	}
	for j := k; j < n; j++ {
		a[j*lda+j] = 1
	}
	for i := k - 1; i >= 0; i-- {
		for i := range work {
			work[i] = 0
		}
		if i < n-1 {
			a[i*lda+i] = 1
			impl.Dlarf(blas.Left, m-i, n-i-1, a[i*lda+i:], lda, tau[i], a[i*lda+i+1:], lda, work)
		}
		if i < m-1 {
			bi.Dscal(m-i-1, -tau[i], a[(i+1)*lda+i:], lda)
		}
		a[i*lda+i] = 1 - tau[i]
		for l := 0; l < i; l++ {
			a[l*lda+i] = 0
		}
	}
}

// Dtrtrs solves a triangular system of the form A * X = B or A^T * X = B. Dtrtrs
// returns whether the solve completed successfully. If A is singular, no solve is performed.
func (impl Implementation) Dtrtrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, a []float64, lda int, b []float64, ldb int) (ok bool) {
	nounit := diag == blas.NonUnit
	if n == 0 {
		return false
	}
	// Check for singularity.
	if nounit {
		for i := 0; i < n; i++ {
			if a[i*lda+i] == 0 {
				return false
			}
		}
	}
	bi := blas64.Implementation()
	bi.Dtrsm(blas.Left, uplo, trans, diag, n, nrhs, 1, a, lda, b, ldb)
	return true
}

// Dorgl2 generates an m×n matrix Q with orthonormal rows defined by the
// first m rows product of elementary reflectors as computed by Dgelqf.
//  Q = H_0 * H_1 * ... * H_{k-1}
// len(tau) >= k, 0 <= k <= m, 0 <= m <= n, len(work) >= m.
// Dorgl2 will panic if these conditions are not met.
//
// Dorgl2 is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dorgl2(m, n, k int, a []float64, lda int, tau, work []float64) {
	checkMatrix(m, n, a, lda)
	if len(tau) < k {
		panic(badTau)
	}
	if k > m {
		panic(kGTM)
	}
	if k > m {
		panic(kGTM)
	}
	if m > n {
		panic(nLTM)
	}
	if len(work) < m {
		panic(badWork)
	}
	if m == 0 {
		return
	}
	bi := blas64.Implementation()
	if k < m {
		for i := k; i < m; i++ {
			for j := 0; j < n; j++ {
				a[i*lda+j] = 0
			}
		}
		for j := k; j < m; j++ {
			a[j*lda+j] = 1
		}
	}
	for i := k - 1; i >= 0; i-- {
		if i < n-1 {
			if i < m-1 {
				a[i*lda+i] = 1
				impl.Dlarf(blas.Right, m-i-1, n-i, a[i*lda+i:], 1, tau[i], a[(i+1)*lda+i:], lda, work)
			}
			bi.Dscal(n-i-1, -tau[i], a[i*lda+i+1:], 1)
		}
		a[i*lda+i] = 1 - tau[i]
		for l := 0; l < i; l++ {
			a[i*lda+l] = 0
		}
	}
}

// Dtrti2 computes the inverse of a triangular matrix, storing the result in place
// into a. This is the BLAS level 2 version of the algorithm.
func (impl Implementation) Dtrti2(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) {
	checkMatrix(n, n, a, lda)
	if uplo != blas.Upper && uplo != blas.Lower {
		panic(badUplo)
	}
	if diag != blas.NonUnit && diag != blas.Unit {
		panic(badDiag)
	}
	bi := blas64.Implementation()

	nonUnit := diag == blas.NonUnit
	// TODO(btracey): Replace this with a row-major ordering.
	if uplo == blas.Upper {
		for j := 0; j < n; j++ {
			var ajj float64
			if nonUnit {
				ajj = 1 / a[j*lda+j]
				a[j*lda+j] = ajj
				ajj *= -1
			} else {
				ajj = -1
			}
			bi.Dtrmv(blas.Upper, blas.NoTrans, diag, j, a, lda, a[j:], lda)
			bi.Dscal(j, ajj, a[j:], lda)
		}
		return
	}
	for j := n - 1; j >= 0; j-- {
		var ajj float64
		if nonUnit {
			ajj = 1 / a[j*lda+j]
			a[j*lda+j] = ajj
			ajj *= -1
		} else {
			ajj = -1
		}
		if j < n-1 {
			bi.Dtrmv(blas.Lower, blas.NoTrans, diag, n-j-1, a[(j+1)*lda+j+1:], lda, a[(j+1)*lda+j:], lda)
			bi.Dscal(n-j-1, ajj, a[(j+1)*lda+j:], lda)
		}
	}
}