// 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)
		}
	}
}

// Dgetrf computes the LU decomposition of the m×n matrix A.
// The LU decomposition is a factorization of A into
//  A = P * L * U
// where P is a permutation matrix, L is a unit lower triangular matrix, and
// U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored
// in place into a.
//
// ipiv is a permutation vector. It indicates that row i of the matrix was
// changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic
// otherwise. ipiv is zero-indexed.
//
// Dgetrf is the blocked version of the algorithm.
//
// Dgetrf returns whether the matrix A is singular. The LU decomposition will
// be computed regardless of the singularity of A, but division by zero
// will occur if the false is returned and the result is used to solve a
// system of equations.
func (impl Implementation) Dgetrf(m, n int, a []float64, lda int, ipiv []int) (ok bool) {
	mn := min(m, n)
	checkMatrix(m, n, a, lda)
	if len(ipiv) < mn {
		panic(badIpiv)
	}
	if m == 0 || n == 0 {
		return
	}
	bi := blas64.Implementation()
	nb := impl.Ilaenv(1, "DGETRF", " ", m, n, -1, -1)
	if nb <= 1 || nb >= min(m, n) {
		// Use the unblocked algorithm.
		return impl.Dgetf2(m, n, a, lda, ipiv)
	}
	ok = true
	for j := 0; j < mn; j += nb {
		jb := min(mn-j, nb)
		blockOk := impl.Dgetf2(m-j, jb, a[j*lda+j:], lda, ipiv[j:])
		if !blockOk {
			ok = false
		}
		for i := j; i <= min(m-1, j+jb-1); i++ {
			ipiv[i] = j + ipiv[i]
		}
		impl.Dlaswp(j, a, lda, j, j+jb-1, ipiv, 1)
		if j+jb < n {
			impl.Dlaswp(n-j-jb, a[j+jb:], lda, j, j+jb-1, ipiv, 1)
			bi.Dtrsm(blas.Left, blas.Lower, blas.NoTrans, blas.Unit,
				jb, n-j-jb, 1,
				a[j*lda+j:], lda,
				a[j*lda+j+jb:], lda)
			if j+jb < m {
				bi.Dgemm(blas.NoTrans, blas.NoTrans, m-j-jb, n-j-jb, jb, -1,
					a[(j+jb)*lda+j:], lda,
					a[j*lda+j+jb:], lda,
					1, a[(j+jb)*lda+j+jb:], lda)
			}
		}
	}
	return ok
}

// Dorg2l generates an m×n matrix Q with orthonormal columns which is defined
// as the last n columns of a product of k elementary reflectors of order m.
//  Q = H_{k-1} * ... * H_1 * H_0
// See Dgelqf for more information. It must be that m >= n >= k.
//
// tau contains the scalar reflectors computed by Dgeqlf. tau must have length
// at least k, and Dorg2l will panic otherwise.
//
// work contains temporary memory, and must have length at least n. Dorg2l will
// panic otherwise.
//
// Dorg2l is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dorg2l(m, n, k int, a []float64, lda int, tau, work []float64) {
	checkMatrix(m, n, a, lda)
	if len(tau) < k {
		panic(badTau)
	}
	if len(work) < n {
		panic(badWork)
	}
	if m < n {
		panic(mLTN)
	}
	if k > n {
		panic(kGTN)
	}
	if n == 0 {
		return
	}

	// Initialize columns 0:n-k to columns of the unit matrix.
	for j := 0; j < n-k; j++ {
		for l := 0; l < m; l++ {
			a[l*lda+j] = 0
		}
		a[(m-n+j)*lda+j] = 1
	}

	bi := blas64.Implementation()
	for i := 0; i < k; i++ {
		ii := n - k + i

		// Apply H_i to A[0:m-k+i, 0:n-k+i] from the left.
		a[(m-n+ii)*lda+ii] = 1
		impl.Dlarf(blas.Left, m-n+ii+1, ii, a[ii:], lda, tau[i], a, lda, work)
		bi.Dscal(m-n+ii, -tau[i], a[ii:], lda)
		a[(m-n+ii)*lda+ii] = 1 - tau[i]

		// Set A[m-k+i:m, n-k+i+1] to zero.
		for l := m - n + ii + 1; l < m; l++ {
			a[l*lda+ii] = 0
		}
	}
}

// Dgetf2 computes the LU decomposition of the m×n matrix A.
// The LU decomposition is a factorization of a into
//  A = P * L * U
// where P is a permutation matrix, L is a unit lower triangular matrix, and
// U is a (usually) non-unit upper triangular matrix. On exit, L and U are stored
// in place into a.
//
// ipiv is a permutation vector. It indicates that row i of the matrix was
// changed with ipiv[i]. ipiv must have length at least min(m,n), and will panic
// otherwise. ipiv is zero-indexed.
//
// Dgetf2 returns whether the matrix A is singular. The LU decomposition will
// be computed regardless of the singularity of A, but division by zero
// will occur if the false is returned and the result is used to solve a
// system of equations.
//
// Dgetf2 is an internal routine. It is exported for testing purposes.
func (Implementation) Dgetf2(m, n int, a []float64, lda int, ipiv []int) (ok bool) {
	mn := min(m, n)
	checkMatrix(m, n, a, lda)
	if len(ipiv) < mn {
		panic(badIpiv)
	}
	if m == 0 || n == 0 {
		return true
	}
	bi := blas64.Implementation()
	sfmin := dlamchS
	ok = true
	for j := 0; j < mn; j++ {
		// Find a pivot and test for singularity.
		jp := j + bi.Idamax(m-j, a[j*lda+j:], lda)
		ipiv[j] = jp
		if a[jp*lda+j] == 0 {
			ok = false
		} else {
			// Swap the rows if necessary.
			if jp != j {
				bi.Dswap(n, a[j*lda:], 1, a[jp*lda:], 1)
			}
			if j < m-1 {
				aj := a[j*lda+j]
				if math.Abs(aj) >= sfmin {
					bi.Dscal(m-j-1, 1/aj, a[(j+1)*lda+j:], lda)
				} else {
					for i := 0; i < m-j-1; i++ {
						a[(j+1)*lda+j] = a[(j+1)*lda+j] / a[lda*j+j]
					}
				}
			}
		}
		if j < mn-1 {
			bi.Dger(m-j-1, n-j-1, -1, a[(j+1)*lda+j:], lda, a[j*lda+j+1:], 1, a[(j+1)*lda+j+1:], lda)
		}
	}
	return ok
}

// Dlarfg generates an elementary reflector for a Householder matrix. It creates
// a real elementary reflector of order n such that
//  H * (alpha) = (beta)
//      (    x)   (   0)
//  H^T * H = I
// H is represented in the form
//  H = 1 - tau * (1; v) * (1 v^T)
// where tau is a real scalar.
//
// On entry, x contains the vector x, on exit it contains v.
func (impl Implementation) Dlarfg(n int, alpha float64, x []float64, incX int) (beta, tau float64) {
	if n < 0 {
		panic(nLT0)
	}
	if n <= 1 {
		return alpha, 0
	}
	checkVector(n-1, x, incX)
	bi := blas64.Implementation()
	xnorm := bi.Dnrm2(n-1, x, incX)
	if xnorm == 0 {
		return alpha, 0
	}
	beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
	safmin := dlamchS / dlamchE
	knt := 0
	if math.Abs(beta) < safmin {
		// xnorm and beta may be innacurate, scale x and recompute.
		rsafmn := 1 / safmin
		for {
			knt++
			bi.Dscal(n-1, rsafmn, x, incX)
			beta *= rsafmn
			alpha *= rsafmn
			if math.Abs(beta) >= safmin {
				break
			}
		}
		xnorm = bi.Dnrm2(n-1, x, incX)
		beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
	}
	tau = (beta - alpha) / beta
	bi.Dscal(n-1, 1/(alpha-beta), x, incX)
	for j := 0; j < knt; j++ {
		beta *= safmin
	}
	return beta, tau
}

func DtrtriTest(t *testing.T, impl Dtrtrier) {
	bi := blas64.Implementation()
	for _, uplo := range []blas.Uplo{blas.Upper} {
		for _, diag := range []blas.Diag{blas.NonUnit, blas.Unit} {
			for _, test := range []struct {
				n, lda int
			}{
				{3, 0},
				{70, 0},
				{200, 0},
				{3, 5},
				{70, 92},
				{200, 205},
			} {
				n := test.n
				lda := test.lda
				if lda == 0 {
					lda = n
				}
				a := make([]float64, n*lda)
				for i := range a {
					a[i] = rand.Float64() + 1 // This keeps the matrices well conditioned.
				}
				aCopy := make([]float64, len(a))
				copy(aCopy, a)
				impl.Dtrtri(uplo, diag, n, a, lda)
				if uplo == blas.Upper {
					for i := 1; i < n; i++ {
						for j := 0; j < i; j++ {
							aCopy[i*lda+j] = 0
							a[i*lda+j] = 0
						}
					}
				} else {
					for i := 1; i < n; i++ {
						for j := i + 1; j < n; j++ {
							aCopy[i*lda+j] = 0
							a[i*lda+j] = 0
						}
					}
				}
				if diag == blas.Unit {
					for i := 0; i < n; i++ {
						a[i*lda+i] = 1
						aCopy[i*lda+i] = 1
					}
				}
				ans := make([]float64, len(a))
				bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, a, lda, aCopy, lda, 0, ans, lda)
				iseye := true
				for i := 0; i < n; i++ {
					for j := 0; j < n; j++ {
						if i == j {
							if math.Abs(ans[i*lda+i]-1) > 1e-4 {
								iseye = false
								break
							}
						} else {
							if math.Abs(ans[i*lda+j]) > 1e-4 {
								iseye = false
								break
							}
						}
					}
				}
				if !iseye {
					t.Errorf("inv(A) * A != I. Upper = %v, unit = %v, n = %v, lda = %v",
						uplo == blas.Upper, diag == blas.Unit, n, lda)
				}
			}
		}
	}
}