Golang blas64.Gemm函数代码示例

func constructH(tau []float64, v blas64.General, store lapack.StoreV, direct lapack.Direct) blas64.General {
	m := v.Rows
	k := v.Cols
	if store == lapack.RowWise {
		m, k = k, m
	}
	h := blas64.General{
		Rows:   m,
		Cols:   m,
		Stride: m,
		Data:   make([]float64, m*m),
	}
	for i := 0; i < m; i++ {
		h.Data[i*m+i] = 1
	}
	for i := 0; i < k; i++ {
		vecData := make([]float64, m)
		if store == lapack.ColumnWise {
			for j := 0; j < m; j++ {
				vecData[j] = v.Data[j*v.Cols+i]
			}
		} else {
			for j := 0; j < m; j++ {
				vecData[j] = v.Data[i*v.Cols+j]
			}
		}
		vec := blas64.Vector{
			Inc:  1,
			Data: vecData,
		}

		hi := blas64.General{
			Rows:   m,
			Cols:   m,
			Stride: m,
			Data:   make([]float64, m*m),
		}
		for i := 0; i < m; i++ {
			hi.Data[i*m+i] = 1
		}
		// hi = I - tau * v * v^T
		blas64.Ger(-tau[i], vec, vec, hi)

		hcopy := blas64.General{
			Rows:   m,
			Cols:   m,
			Stride: m,
			Data:   make([]float64, m*m),
		}
		copy(hcopy.Data, h.Data)
		if direct == lapack.Forward {
			// H = H * H_I in forward mode
			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hcopy, hi, 0, h)
		} else {
			// H = H_I * H in backward mode
			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hi, hcopy, 0, h)
		}
	}
	return h
}

func (c *ConvLayer) convolveR(v autofunc.RVector, in, inR linalg.Vector, out *Tensor3) {
	inMat := c.inputToMatrix(in)
	inMatR := c.inputToMatrix(inR)
	filterMat := blas64.General{
		Rows:   c.FilterCount,
		Cols:   inMat.Cols,
		Stride: inMat.Stride,
		Data:   c.FilterVar.Vector,
	}
	outMat := blas64.General{
		Rows:   out.Width * out.Height,
		Cols:   out.Depth,
		Stride: out.Depth,
		Data:   out.Data,
	}
	blas64.Gemm(blas.NoTrans, blas.Trans, 1, inMatR, filterMat, 0, outMat)
	if filterRV, ok := v[c.FilterVar]; ok {
		filterMatR := blas64.General{
			Rows:   c.FilterCount,
			Cols:   inMat.Cols,
			Stride: inMat.Stride,
			Data:   filterRV,
		}
		blas64.Gemm(blas.NoTrans, blas.Trans, 1, inMat, filterMatR, 1, outMat)
	}

	if biasRV, ok := v[c.Biases]; ok {
		biasVec := blas64.Vector{Inc: 1, Data: biasRV}
		for i := 0; i < len(out.Data); i += outMat.Cols {
			outRow := out.Data[i : i+outMat.Cols]
			outVec := blas64.Vector{Inc: 1, Data: outRow}
			blas64.Axpy(len(outRow), 1, biasVec, outVec)
		}
	}
}

func (c *convLayerResult) propagateSingle(input, upstream, downstream linalg.Vector,
	grad autofunc.Gradient) {
	upstreamMat := blas64.General{
		Rows:   c.Layer.OutputWidth() * c.Layer.OutputHeight(),
		Cols:   c.Layer.OutputDepth(),
		Stride: c.Layer.OutputDepth(),
		Data:   upstream,
	}

	if downstream != nil {
		inDeriv := c.Layer.inputToMatrix(input)
		filterMat := blas64.General{
			Rows:   len(c.Layer.Filters),
			Cols:   c.Layer.FilterWidth * c.Layer.FilterHeight * c.Layer.InputDepth,
			Stride: c.Layer.FilterWidth * c.Layer.FilterHeight * c.Layer.InputDepth,
			Data:   c.Layer.FilterVar.Vector,
		}
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, upstreamMat, filterMat, 0, inDeriv)
		flattened := NewTensor3Col(c.Layer.InputWidth, c.Layer.InputHeight,
			c.Layer.InputDepth, inDeriv.Data, c.Layer.FilterWidth,
			c.Layer.FilterHeight, c.Layer.Stride)
		copy(downstream, flattened.Data)
	}

	if filterGrad, ok := grad[c.Layer.FilterVar]; ok {
		inMatrix := c.Layer.inputToMatrix(input)
		destMat := blas64.General{
			Rows:   len(c.Layer.Filters),
			Cols:   c.Layer.FilterWidth * c.Layer.FilterHeight * c.Layer.InputDepth,
			Stride: c.Layer.FilterWidth * c.Layer.FilterHeight * c.Layer.InputDepth,
			Data:   filterGrad,
		}
		blas64.Gemm(blas.Trans, blas.NoTrans, 1, upstreamMat, inMatrix, 1, destMat)
	}
}

func testDlarfx(t *testing.T, impl Dlarfxer, side blas.Side, m, n, extra int, rnd *rand.Rand) {
	const tol = 1e-13

	c := randomGeneral(m, n, n+extra, rnd)
	cWant := randomGeneral(m, n, n+extra, rnd)
	tau := rnd.NormFloat64()

	var (
		v []float64
		h blas64.General
	)
	if side == blas.Left {
		v = randomSlice(m, rnd)
		h = eye(m, m+extra)
	} else {
		v = randomSlice(n, rnd)
		h = eye(n, n+extra)
	}
	blas64.Ger(-tau, blas64.Vector{Inc: 1, Data: v}, blas64.Vector{Inc: 1, Data: v}, h)
	if side == blas.Left {
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, c, 0, cWant)
	} else {
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, c, h, 0, cWant)
	}

	var work []float64
	if h.Rows > 10 {
		// Allocate work only if H has order > 10.
		if side == blas.Left {
			work = make([]float64, n)
		} else {
			work = make([]float64, m)
		}
	}

	impl.Dlarfx(side, m, n, v, tau, c.Data, c.Stride, work)

	prefix := fmt.Sprintf("Case side=%v, m=%v, n=%v, extra=%v", side, m, n, extra)

	// Check any invalid modifications of c.
	if !generalOutsideAllNaN(c) {
		t.Errorf("%v: out-of-range write to C\n%v", prefix, c.Data)
	}

	if !equalApproxGeneral(c, cWant, tol) {
		t.Errorf("%v: unexpected C\n%v", prefix, c.Data)
	}
}

// Mul takes the matrix product of a and b, placing the result in the receiver.
//
// See the Muler interface for more information.
func (m *Dense) Mul(a, b Matrix) {
	ar, ac := a.Dims()
	br, bc := b.Dims()

	if ac != br {
		panic(ErrShape)
	}

	m.reuseAs(ar, bc)
	var w *Dense
	if m != a && m != b {
		w = m
	} else {
		w = getWorkspace(ar, bc, false)
		defer func() {
			m.Copy(w)
			putWorkspace(w)
		}()
	}

	if a, ok := a.(RawMatrixer); ok {
		if b, ok := b.(RawMatrixer); ok {
			amat, bmat := a.RawMatrix(), b.RawMatrix()
			blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, amat, bmat, 0, w.Mat)
			return
		}
	}

	if a, ok := a.(Vectorer); ok {
		if b, ok := b.(Vectorer); ok {
			row := make([]float64, ac)
			col := make([]float64, br)
			for r := 0; r < ar; r++ {
				dataTmp := w.Mat.Data[r*w.Mat.Stride : r*w.Mat.Stride+bc]
				for c := 0; c < bc; c++ {
					dataTmp[c] = blas64.Dot(ac,
						blas64.Vector{Inc: 1, Data: a.Row(row, r)},
						blas64.Vector{Inc: 1, Data: b.Col(col, c)},
					)
				}
			}
			return
		}
	}

	row := make([]float64, ac)
	for r := 0; r < ar; r++ {
		for i := range row {
			row[i] = a.At(r, i)
		}
		for c := 0; c < bc; c++ {
			var v float64
			for i, e := range row {
				v += e * b.At(i, c)
			}
			w.Mat.Data[r*w.Mat.Stride+c] = v
		}
	}
}

func testDorghr(t *testing.T, impl Dorghrer, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
	const tol = 1e-14

	// Construct the matrix A with elementary reflectors and scalar factors tau.
	a := randomGeneral(n, n, n+extra, rnd)
	var tau []float64
	if n > 1 {
		tau = nanSlice(n - 1)
	}
	work := nanSlice(max(1, n)) // Minimum work for Dgehrd.
	impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))

	// Extract Q for later comparison.
	q := eye(n, n)
	qCopy := cloneGeneral(q)
	for j := ilo; j < ihi; j++ {
		h := eye(n, n)
		v := blas64.Vector{
			Inc:  1,
			Data: make([]float64, n),
		}
		v.Data[j+1] = 1
		for i := j + 2; i < ihi+1; i++ {
			v.Data[i] = a.Data[i*a.Stride+j]
		}
		blas64.Ger(-tau[j], v, v, h)
		copy(qCopy.Data, q.Data)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
	}

	if optwork {
		work = nanSlice(1)
		impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, -1)
		work = nanSlice(int(work[0]))
	} else {
		work = nanSlice(max(1, ihi-ilo))
	}
	impl.Dorghr(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))

	prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v, optwork=%v", n, ilo, ihi, extra, optwork)
	if !generalOutsideAllNaN(a) {
		t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
	}
	if !isOrthonormal(a) {
		t.Errorf("%v: A is not orthogonal\n%v", prefix, a.Data)
	}
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			aij := a.Data[i*a.Stride+j]
			qij := q.Data[i*q.Stride+j]
			if math.Abs(aij-qij) > tol {
				t.Errorf("%v: unexpected value of A[%v,%v]. want %v, got %v", prefix, i, j, qij, aij)
			}
		}
	}
}

// QFromQR extracts the m×m orthonormal matrix Q from a QR decomposition.
func (m *Dense) QFromQR(qr *QR) {
	r, c := qr.qr.Dims()
	m.reuseAs(r, r)

	// Set Q = I.
	for i := 0; i < r; i++ {
		for j := 0; j < i; j++ {
			m.mat.Data[i*m.mat.Stride+j] = 0
		}
		m.mat.Data[i*m.mat.Stride+i] = 1
		for j := i + 1; j < r; j++ {
			m.mat.Data[i*m.mat.Stride+j] = 0
		}
	}

	// Construct Q from the elementary reflectors.
	h := blas64.General{
		Rows:   r,
		Cols:   r,
		Stride: r,
		Data:   make([]float64, r*r),
	}
	qCopy := getWorkspace(r, r, false)
	v := blas64.Vector{
		Inc:  1,
		Data: make([]float64, r),
	}
	for i := 0; i < c; i++ {
		// Set h = I.
		for i := range h.Data {
			h.Data[i] = 0
		}
		for j := 0; j < r; j++ {
			h.Data[j*r+j] = 1
		}

		// Set the vector data as the elementary reflector.
		for j := 0; j < i; j++ {
			v.Data[j] = 0
		}
		v.Data[i] = 1
		for j := i + 1; j < r; j++ {
			v.Data[j] = qr.qr.mat.Data[j*qr.qr.mat.Stride+i]
		}

		// Compute the multiplication matrix.
		blas64.Ger(-qr.tau[i], v, v, h)
		qCopy.Copy(m)
		blas64.Gemm(blas.NoTrans, blas.NoTrans,
			1, qCopy.mat, h,
			0, m.mat)
	}
}

// QFromLQ extracts the n×n orthonormal matrix Q from an LQ decomposition.
func (m *Dense) QFromLQ(lq *LQ) {
	r, c := lq.lq.Dims()
	m.reuseAs(c, c)

	// Set Q = I.
	for i := 0; i < c; i++ {
		for j := 0; j < i; j++ {
			m.mat.Data[i*m.mat.Stride+j] = 0
		}
		m.mat.Data[i*m.mat.Stride+i] = 1
		for j := i + 1; j < c; j++ {
			m.mat.Data[i*m.mat.Stride+j] = 0
		}
	}

	// Construct Q from the elementary reflectors.
	h := blas64.General{
		Rows:   c,
		Cols:   c,
		Stride: c,
		Data:   make([]float64, c*c),
	}
	qCopy := getWorkspace(c, c, false)
	v := blas64.Vector{
		Inc:  1,
		Data: make([]float64, c),
	}
	for i := 0; i < r; i++ {
		// Set h = I.
		for i := range h.Data {
			h.Data[i] = 0
		}
		for j := 0; j < c; j++ {
			h.Data[j*c+j] = 1
		}

		// Set the vector data as the elementary reflector.
		for j := 0; j < i; j++ {
			v.Data[j] = 0
		}
		v.Data[i] = 1
		for j := i + 1; j < c; j++ {
			v.Data[j] = lq.lq.mat.Data[i*lq.lq.mat.Stride+j]
		}

		// Compute the multiplication matrix.
		blas64.Ger(-lq.tau[i], v, v, h)
		qCopy.Copy(m)
		blas64.Gemm(blas.NoTrans, blas.NoTrans,
			1, h, qCopy.mat,
			0, m.mat)
	}
}

// Eval returns a matrix literal.
func (m1 *Mul) Eval() MatrixLiteral {

	// This should be replaced with a call to Eval on each side, and then a type
	// switch to handle the various matrix literals.

	lm := m1.Left.Eval()
	rm := m1.Right.Eval()

	left := lm.AsGeneral()
	right := rm.AsGeneral()
	r, c := m1.Dims()
	m := blas64.General{
		Rows:   r,
		Cols:   c,
		Stride: c,
		Data:   make([]float64, r*c),
	}
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, left, right, 0, m)
	return &General{m}
}

func (c *ConvLayer) convolve(in linalg.Vector, out *Tensor3) {
	inMat := c.inputToMatrix(in)
	filterMat := blas64.General{
		Rows:   c.FilterCount,
		Cols:   inMat.Cols,
		Stride: inMat.Stride,
		Data:   c.FilterVar.Vector,
	}
	outMat := blas64.General{
		Rows:   out.Width * out.Height,
		Cols:   out.Depth,
		Stride: out.Depth,
		Data:   out.Data,
	}
	blas64.Gemm(blas.NoTrans, blas.Trans, 1, inMat, filterMat, 0, outMat)

	biasVec := blas64.Vector{Inc: 1, Data: c.Biases.Vector}
	for i := 0; i < len(out.Data); i += outMat.Cols {
		outRow := out.Data[i : i+outMat.Cols]
		outVec := blas64.Vector{Inc: 1, Data: outRow}
		blas64.Axpy(len(outRow), 1, biasVec, outVec)
	}
}

// dlatrdCheckDecomposition checks that the first nb rows have been successfully
// reduced.
func dlatrdCheckDecomposition(t *testing.T, uplo blas.Uplo, n, nb int, e, tau, a []float64, lda int, aGen, q blas64.General) bool {
	// Compute Q^T * A * Q.
	tmp := blas64.General{
		Rows:   n,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, n*n),
	}

	ans := blas64.General{
		Rows:   n,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, n*n),
	}

	blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aGen, 0, tmp)
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, q, 0, ans)

	// Compare with T.
	if uplo == blas.Upper {
		for i := n - 1; i >= n-nb; i-- {
			for j := 0; j < n; j++ {
				v := ans.Data[i*ans.Stride+j]
				switch {
				case i == j:
					if math.Abs(v-a[i*lda+j]) > 1e-10 {
						return false
					}
				case i == j-1:
					if math.Abs(a[i*lda+j]-1) > 1e-10 {
						return false
					}
					if math.Abs(v-e[i]) > 1e-10 {
						return false
					}
				case i == j+1:
				default:
					if math.Abs(v) > 1e-10 {
						return false
					}
				}
			}
		}
	} else {
		for i := 0; i < nb; i++ {
			for j := 0; j < n; j++ {
				v := ans.Data[i*ans.Stride+j]
				switch {
				case i == j:
					if math.Abs(v-a[i*lda+j]) > 1e-10 {
						return false
					}
				case i == j-1:
				case i == j+1:
					if math.Abs(a[i*lda+j]-1) > 1e-10 {
						return false
					}
					if math.Abs(v-e[i-1]) > 1e-10 {
						return false
					}
				default:
					if math.Abs(v) > 1e-10 {
						return false
					}
				}
			}
		}
	}
	return true
}

//.........这里部分代码省略.........
						a := make([]float64, ma*lda)
						for i := 0; i < ma; i++ {
							for j := 0; j < lda; j++ {
								a[i*lda+j] = rnd.Float64()
							}
						}
						k := min(ma, na)

						// H is always ma x ma
						var m, n, rowsWork int
						switch {
						default:
							panic("not implemented")
						case side == blas.Left:
							m = test.ma
							n = test.cdim
							rowsWork = n
						case side == blas.Right:
							m = test.cdim
							n = test.ma
							rowsWork = m
						}

						// Use dgeqr2 to find the v vectors
						tau := make([]float64, na)
						work := make([]float64, na)
						impl.Dgeqr2(ma, k, a, lda, tau, work)

						// Correct the v vectors based on the direct and store
						vMatTmp := extractVMat(ma, na, a, lda, lapack.Forward, lapack.ColumnWise)
						vMat := constructVMat(vMatTmp, store, direct)
						v := vMat.Data
						ldv := vMat.Stride

						// Use dlarft to find the t vector
						ldt := test.ldt
						if ldt == 0 {
							ldt = k
						}
						tm := make([]float64, k*ldt)

						impl.Dlarft(direct, store, ma, k, v, ldv, tau, tm, ldt)

						// Generate c matrix
						ldc := test.ldc
						if ldc == 0 {
							ldc = n
						}
						c := make([]float64, m*ldc)
						for i := 0; i < m; i++ {
							for j := 0; j < ldc; j++ {
								c[i*ldc+j] = rnd.Float64()
							}
						}
						cCopy := make([]float64, len(c))
						copy(cCopy, c)

						ldwork := k
						work = make([]float64, rowsWork*k)

						// Call Dlarfb with this information
						impl.Dlarfb(side, trans, direct, store, m, n, k, v, ldv, tm, ldt, c, ldc, work, ldwork)

						h := constructH(tau, vMat, store, direct)

						cMat := blas64.General{
							Rows:   m,
							Cols:   n,
							Stride: ldc,
							Data:   make([]float64, m*ldc),
						}
						copy(cMat.Data, cCopy)
						ans := blas64.General{
							Rows:   m,
							Cols:   n,
							Stride: ldc,
							Data:   make([]float64, m*ldc),
						}
						copy(ans.Data, cMat.Data)
						switch {
						default:
							panic("not implemented")
						case side == blas.Left && trans == blas.NoTrans:
							blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, cMat, 0, ans)
						case side == blas.Left && trans == blas.Trans:
							blas64.Gemm(blas.Trans, blas.NoTrans, 1, h, cMat, 0, ans)
						case side == blas.Right && trans == blas.NoTrans:
							blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMat, h, 0, ans)
						case side == blas.Right && trans == blas.Trans:
							blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMat, h, 0, ans)
						}
						if !floats.EqualApprox(ans.Data, c, 1e-14) {
							t.Errorf("Cas %v mismatch. Want %v, got %v.", cas, ans.Data, c)
						}
					}
				}
			}
		}
	}
}

//.........这里部分代码省略.........

			lastr: 0,
			lastc: 1,

			tau: 2,
		},
		{
			m:   10,
			n:   10,
			ldc: 10,

			incv:  4,
			lastv: 6,

			lastr: 9,
			lastc: 8,

			tau: 2,
		},
	} {
		// Construct a random matrix.
		c := make([]float64, test.ldc*test.m)
		for i := 0; i <= test.lastr; i++ {
			for j := 0; j <= test.lastc; j++ {
				c[i*test.ldc+j] = rand.Float64()
			}
		}
		cCopy := make([]float64, len(c))
		copy(cCopy, c)
		cCopy2 := make([]float64, len(c))
		copy(cCopy2, c)

		// Test with side right.
		sz := max(test.m, test.n) // so v works for both right and left side.
		v := make([]float64, test.incv*sz+1)
		// Fill with nonzero entries up until lastv.
		for i := 0; i <= test.lastv; i++ {
			v[i*test.incv] = rand.Float64()
		}
		// Construct h explicitly to compare.
		h := make([]float64, test.n*test.n)
		for i := 0; i < test.n; i++ {
			h[i*test.n+i] = 1
		}
		hMat := blas64.General{
			Rows:   test.n,
			Cols:   test.n,
			Stride: test.n,
			Data:   h,
		}
		vVec := blas64.Vector{
			Inc:  test.incv,
			Data: v,
		}
		blas64.Ger(-test.tau, vVec, vVec, hMat)

		// Apply multiplication (2nd copy is to avoid aliasing).
		cMat := blas64.General{
			Rows:   test.m,
			Cols:   test.n,
			Stride: test.ldc,
			Data:   cCopy,
		}
		cMat2 := blas64.General{
			Rows:   test.m,
			Cols:   test.n,
			Stride: test.ldc,
			Data:   cCopy2,
		}
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMat2, hMat, 0, cMat)

		// cMat now stores the true answer. Compare with the function call.
		work := make([]float64, sz)
		impl.Dlarf(blas.Right, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
		if !floats.EqualApprox(c, cMat.Data, 1e-14) {
			t.Errorf("Dlarf mismatch right, case %v. Want %v, got %v", i, cMat.Data, c)
		}

		// Test on the left side.
		copy(c, cCopy2)
		copy(cCopy, c)
		// Construct h.
		h = make([]float64, test.m*test.m)
		for i := 0; i < test.m; i++ {
			h[i*test.m+i] = 1
		}
		hMat = blas64.General{
			Rows:   test.m,
			Cols:   test.m,
			Stride: test.m,
			Data:   h,
		}
		blas64.Ger(-test.tau, vVec, vVec, hMat)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hMat, cMat2, 0, cMat)
		impl.Dlarf(blas.Left, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
		if !floats.EqualApprox(c, cMat.Data, 1e-14) {
			t.Errorf("Dlarf mismatch left, case %v. Want %v, got %v", i, cMat.Data, c)
		}
	}
}

// checkPLU checks that the PLU factorization contained in factorize matches
// the original matrix contained in original.
func checkPLU(t *testing.T, ok bool, m, n, lda int, ipiv []int, factorized, original []float64, tol float64, print bool) {
	var hasZeroDiagonal bool
	for i := 0; i < min(m, n); i++ {
		if factorized[i*lda+i] == 0 {
			hasZeroDiagonal = true
			break
		}
	}
	if hasZeroDiagonal && ok {
		t.Error("Has a zero diagonal but returned ok")
	}
	if !hasZeroDiagonal && !ok {
		t.Error("Non-zero diagonal but returned !ok")
	}

	// Check that the LU decomposition is correct.
	mn := min(m, n)
	l := make([]float64, m*mn)
	ldl := mn
	u := make([]float64, mn*n)
	ldu := n
	for i := 0; i < m; i++ {
		for j := 0; j < n; j++ {
			v := factorized[i*lda+j]
			switch {
			case i == j:
				l[i*ldl+i] = 1
				u[i*ldu+i] = v
			case i > j:
				l[i*ldl+j] = v
			case i < j:
				u[i*ldu+j] = v
			}
		}
	}

	LU := blas64.General{
		Rows:   m,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, m*n),
	}
	U := blas64.General{
		Rows:   mn,
		Cols:   n,
		Stride: ldu,
		Data:   u,
	}
	L := blas64.General{
		Rows:   m,
		Cols:   mn,
		Stride: ldl,
		Data:   l,
	}
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, L, U, 0, LU)

	p := make([]float64, m*m)
	ldp := m
	for i := 0; i < m; i++ {
		p[i*ldp+i] = 1
	}
	for i := len(ipiv) - 1; i >= 0; i-- {
		v := ipiv[i]
		blas64.Swap(m, blas64.Vector{1, p[i*ldp:]}, blas64.Vector{1, p[v*ldp:]})
	}
	P := blas64.General{
		Rows:   m,
		Cols:   m,
		Stride: m,
		Data:   p,
	}
	aComp := blas64.General{
		Rows:   m,
		Cols:   n,
		Stride: lda,
		Data:   make([]float64, m*lda),
	}
	copy(aComp.Data, factorized)
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, P, LU, 0, aComp)
	if !floats.EqualApprox(aComp.Data, original, tol) {
		if print {
			t.Errorf("PLU multiplication does not match original matrix.\nWant: %v\nGot: %v", original, aComp.Data)
			return
		}
		t.Error("PLU multiplication does not match original matrix.")
	}
}

func Dgelq2Test(t *testing.T, impl Dgelq2er) {
	for c, test := range []struct {
		m, n, lda int
	}{
		{1, 1, 0},
		{2, 2, 0},
		{3, 2, 0},
		{2, 3, 0},
		{1, 12, 0},
		{2, 6, 0},
		{3, 4, 0},
		{4, 3, 0},
		{6, 2, 0},
		{1, 12, 0},
		{1, 1, 20},
		{2, 2, 20},
		{3, 2, 20},
		{2, 3, 20},
		{1, 12, 20},
		{2, 6, 20},
		{3, 4, 20},
		{4, 3, 20},
		{6, 2, 20},
		{1, 12, 20},
	} {
		n := test.n
		m := test.m
		lda := test.lda
		if lda == 0 {
			lda = test.n
		}
		k := min(m, n)
		tau := make([]float64, k)
		for i := range tau {
			tau[i] = rand.Float64()
		}
		work := make([]float64, m)
		for i := range work {
			work[i] = rand.Float64()
		}
		a := make([]float64, m*lda)
		for i := 0; i < m*lda; i++ {
			a[i] = rand.Float64()
		}
		aCopy := make([]float64, len(a))
		copy(aCopy, a)
		impl.Dgelq2(m, n, a, lda, tau, work)

		Q := constructQ("LQ", m, n, a, lda, tau)

		// Check that Q is orthonormal
		for i := 0; i < Q.Rows; i++ {
			nrm := blas64.Nrm2(Q.Cols, blas64.Vector{Inc: 1, Data: Q.Data[i*Q.Stride:]})
			if math.Abs(nrm-1) > 1e-14 {
				t.Errorf("Q not normal. Norm is %v", nrm)
			}
			for j := 0; j < i; j++ {
				dot := blas64.Dot(Q.Rows,
					blas64.Vector{Inc: 1, Data: Q.Data[i*Q.Stride:]},
					blas64.Vector{Inc: 1, Data: Q.Data[j*Q.Stride:]},
				)
				if math.Abs(dot) > 1e-14 {
					t.Errorf("Q not orthogonal. Dot is %v", dot)
				}
			}
		}

		L := blas64.General{
			Rows:   m,
			Cols:   n,
			Stride: n,
			Data:   make([]float64, m*n),
		}
		for i := 0; i < m; i++ {
			for j := 0; j <= min(i, n-1); j++ {
				L.Data[i*L.Stride+j] = a[i*lda+j]
			}
		}

		ans := blas64.General{
			Rows:   m,
			Cols:   n,
			Stride: lda,
			Data:   make([]float64, m*lda),
		}
		copy(ans.Data, aCopy)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, L, Q, 0, ans)
		if !floats.EqualApprox(aCopy, ans.Data, 1e-14) {
			t.Errorf("Case %v, LQ mismatch. Want %v, got %v.", c, aCopy, ans.Data)
		}
	}
}