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

// Dot returns the sum of the element-wise product of a and b.
// Dot panics if the matrix sizes are unequal.
func Dot(a, b *Vector) float64 {
	la := a.Len()
	lb := b.Len()
	if la != lb {
		panic(matrix.ErrShape)
	}
	return blas64.Dot(la, a.mat, b.mat)
}

// Dot returns the sum of the element-wise products of the elements of a and b.
// Dot panics if the matrix sizes are unequal.
func Dot(a, b Matrix) float64 {
	r, c := a.Dims()
	rb, cb := b.Dims()
	if r != rb || c != cb {
		panic(matrix.ErrShape)
	}
	aU, aTrans := untranspose(a)
	bU, bTrans := untranspose(b)
	if rma, ok := aU.(RawVectorer); ok {
		if rmb, ok := bU.(RawVectorer); ok {
			if c > r {
				r = c
			}
			return blas64.Dot(r, rma.RawVector(), rmb.RawVector())
		}
	}
	var sum float64
	if rma, ok := aU.(RawMatrixer); ok {
		if rmb, ok := bU.(RawMatrixer); ok {
			ra := rma.RawMatrix()
			rb := rmb.RawMatrix()
			if aTrans == bTrans {
				for i := 0; i < ra.Rows; i++ {
					sum += blas64.Dot(ra.Cols,
						blas64.Vector{Inc: 1, Data: ra.Data[i*ra.Stride:]},
						blas64.Vector{Inc: 1, Data: rb.Data[i*rb.Stride:]},
					)
				}
				return sum
			}
			for i := 0; i < ra.Rows; i++ {
				sum += blas64.Dot(ra.Cols,
					blas64.Vector{Inc: 1, Data: ra.Data[i*ra.Stride:]},
					blas64.Vector{Inc: rb.Stride, Data: rb.Data[i:]},
				)
			}
			return sum
		}
	}
	for i := 0; i < r; i++ {
		for j := 0; j < c; j++ {
			sum += a.At(i, j) * b.At(i, j)
		}
	}
	return sum
}

func newSimilarityCircuit(uVal, uGrad, vVal, vGrad []float64) *similarityCircuit {
	s := similarityCircuit{
		UVal:  uVal,
		UGrad: uGrad,
		VVal:  vVal,
		VGrad: vGrad,
	}
	u := blas64.Vector{Inc: 1, Data: uVal}
	v := blas64.Vector{Inc: 1, Data: vVal}
	s.UV = blas64.Dot(len(uVal), u, v)
	s.Unorm = blas64.Nrm2(len(uVal), u)
	s.Vnorm = blas64.Nrm2(len(vVal), v)
	s.TopVal = s.UV / (s.Unorm * s.Vnorm)
	return &s
}

// hasOrthornormalColumns checks that the columns of a are orthonormal.
func hasOrthonormalColumns(m, n int, a []float64, lda int) bool {
	for i := 0; i < n; i++ {
		for j := i; j < n; j++ {
			dot := blas64.Dot(m,
				blas64.Vector{Inc: lda, Data: a[i:]},
				blas64.Vector{Inc: lda, Data: a[j:]},
			)
			if i == j {
				if math.Abs(dot-1) > 1e-10 {
					return false
				}
			} else {
				if math.Abs(dot) > 1e-10 {
					return false
				}
			}
		}
	}
	return true
}

// isOrthonormal checks that a general matrix is orthonormal.
func isOrthonormal(q blas64.General) bool {
	n := q.Rows
	for i := 0; i < n; i++ {
		for j := i; j < n; j++ {
			dot := blas64.Dot(n,
				blas64.Vector{Inc: 1, Data: q.Data[i*q.Stride:]},
				blas64.Vector{Inc: 1, Data: q.Data[j*q.Stride:]},
			)
			if i == j {
				if math.Abs(dot-1) > 1e-10 {
					return false
				}
			} else {
				if math.Abs(dot) > 1e-10 {
					return false
				}
			}
		}
	}
	return true
}

func isOrthonormal(q *Dense, tol float64) bool {
	m, n := q.Dims()
	if m != n {
		return false
	}
	for i := 0; i < m; i++ {
		for j := i; j < m; j++ {
			dot := blas64.Dot(m,
				blas64.Vector{Inc: 1, Data: q.mat.Data[i*q.mat.Stride:]},
				blas64.Vector{Inc: 1, Data: q.mat.Data[j*q.mat.Stride:]},
			)
			// Dot product should be 1 if i == j and 0 otherwise.
			if i == j && math.Abs(dot-1) > tol {
				return false
			}
			if i != j && math.Abs(dot) > tol {
				return false
			}
		}
	}
	return true
}

// MulVec computes a * b if trans == false and a^T * b if trans == true. The
// result is stored into the reciever. MulVec panics if the number of columns in
// a does not equal the number of rows in b.
func (m *Vector) MulVec(a Matrix, trans bool, b *Vector) {
	ar, ac := a.Dims()
	br, _ := b.Dims()
	if trans {
		if ar != br {
			panic(ErrShape)
		}
	} else {
		if ac != br {
			panic(ErrShape)
		}
	}

	var w Vector
	if m != a && m != b {
		w = *m
	}
	if w.n == 0 {
		if trans {
			w.mat.Data = use(w.mat.Data, ac)
		} else {
			w.mat.Data = use(w.mat.Data, ar)
		}

		w.mat.Inc = 1
		w.n = ar
	} else {
		if trans {
			if ac != w.n {
				panic(ErrShape)
			}
		} else {
			if ar != w.n {
				panic(ErrShape)
			}
		}
	}

	switch a := a.(type) {
	case RawSymmetricer:
		amat := a.RawSymmetric()
		blas64.Symv(1, amat, b.mat, 0, w.mat)
		*m = w
		return
	case RawMatrixer:
		amat := a.RawMatrix()
		t := blas.NoTrans
		if trans {
			t = blas.Trans
		}
		blas64.Gemv(t, 1, amat, b.mat, 0, w.mat)
		*m = w
		return
	case Vectorer:
		row := make([]float64, ac)
		for r := 0; r < ar; r++ {
			w.mat.Data[r*m.mat.Inc] = blas64.Dot(ac,
				blas64.Vector{Inc: 1, Data: a.Row(row, r)},
				b.mat,
			)
		}
		*m = w
		return
	default:
		row := make([]float64, ac)
		for r := 0; r < ar; r++ {
			for i := range row {
				row[i] = a.At(r, i)
			}
			var v float64
			for i, e := range row {
				v += e * b.mat.Data[i*b.mat.Inc]
			}
			w.mat.Data[r*m.mat.Inc] = v
		}
		*m = w
		return
	}
}

// MulVec computes a * b. The result is stored into the receiver.
// MulVec panics if the number of columns in a does not equal the number of rows in b.
func (v *Vector) MulVec(a Matrix, b *Vector) {
	r, c := a.Dims()
	br := b.Len()
	if c != br {
		panic(ErrShape)
	}
	a, trans := untranspose(a)
	ar, ac := a.Dims()
	v.reuseAs(r)
	var restore func()
	if v == a {
		v, restore = v.isolatedWorkspace(a.(*Vector))
		defer restore()
	} else if v == b {
		v, restore = v.isolatedWorkspace(b)
		defer restore()
	}

	switch a := a.(type) {
	case *Vector:
		if a.Len() == 1 {
			// {1,1} x {1,n}
			av := a.At(0, 0)
			for i := 0; i < b.Len(); i++ {
				v.mat.Data[i*v.mat.Inc] = av * b.mat.Data[i*b.mat.Inc]
			}
			return
		}
		if b.Len() == 1 {
			// {1,n} x {1,1}
			bv := b.At(0, 0)
			for i := 0; i < a.Len(); i++ {
				v.mat.Data[i*v.mat.Inc] = bv * a.mat.Data[i*a.mat.Inc]
			}
			return
		}
		// {n,1} x {1,n}
		var sum float64
		for i := 0; i < c; i++ {
			sum += a.At(i, 0) * b.At(i, 0)
		}
		v.SetVec(0, sum)
		return
	case RawSymmetricer:
		amat := a.RawSymmetric()
		blas64.Symv(1, amat, b.mat, 0, v.mat)
	case RawTriangular:
		v.CopyVec(b)
		amat := a.RawTriangular()
		ta := blas.NoTrans
		if trans {
			ta = blas.Trans
		}
		blas64.Trmv(ta, amat, v.mat)
	case RawMatrixer:
		amat := a.RawMatrix()
		t := blas.NoTrans
		if trans {
			t = blas.Trans
		}
		blas64.Gemv(t, 1, amat, b.mat, 0, v.mat)
	case Vectorer:
		if trans {
			col := make([]float64, ar)
			for c := 0; c < ac; c++ {
				v.mat.Data[c*v.mat.Inc] = blas64.Dot(ar,
					blas64.Vector{Inc: 1, Data: a.Col(col, c)},
					b.mat,
				)
			}
		} else {
			row := make([]float64, ac)
			for r := 0; r < ar; r++ {
				v.mat.Data[r*v.mat.Inc] = blas64.Dot(ac,
					blas64.Vector{Inc: 1, Data: a.Row(row, r)},
					b.mat,
				)
			}
		}
	default:
		if trans {
			col := make([]float64, ar)
			for c := 0; c < ac; c++ {
				for i := range col {
					col[i] = a.At(i, c)
				}
				var f float64
				for i, e := range col {
					f += e * b.mat.Data[i*b.mat.Inc]
				}
				v.mat.Data[c*v.mat.Inc] = f
			}
		} else {
			row := make([]float64, ac)
			for r := 0; r < ar; r++ {
				for i := range row {
					row[i] = a.At(r, i)
				}
//.........这里部分代码省略.........

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

// MulVec computes a * b if trans == false and a^T * b if trans == true. The
// result is stored into the receiver. MulVec panics if the number of columns in
// a does not equal the number of rows in b.
func (v *Vector) MulVec(a Matrix, trans bool, b *Vector) {
	ar, ac := a.Dims()
	br := b.Len()
	if trans {
		if ar != br {
			panic(ErrShape)
		}
	} else {
		if ac != br {
			panic(ErrShape)
		}
	}

	var w Vector
	if v != a && v != b {
		w = *v
	}
	if w.n == 0 {
		if trans {
			w.mat.Data = use(w.mat.Data, ac)
		} else {
			w.mat.Data = use(w.mat.Data, ar)
		}

		w.mat.Inc = 1
		w.n = ar
		if trans {
			w.n = ac
		}
	} else {
		if trans {
			if ac != w.n {
				panic(ErrShape)
			}
		} else {
			if ar != w.n {
				panic(ErrShape)
			}
		}
	}

	switch a := a.(type) {
	case RawSymmetricer:
		amat := a.RawSymmetric()
		blas64.Symv(1, amat, b.mat, 0, w.mat)
	case RawTriangular:
		w.CopyVec(b)
		amat := a.RawTriangular()
		ta := blas.NoTrans
		if trans {
			ta = blas.Trans
		}
		blas64.Trmv(ta, amat, w.mat)
	case RawMatrixer:
		amat := a.RawMatrix()
		t := blas.NoTrans
		if trans {
			t = blas.Trans
		}
		blas64.Gemv(t, 1, amat, b.mat, 0, w.mat)
	case Vectorer:
		if trans {
			col := make([]float64, ar)
			for c := 0; c < ac; c++ {
				w.mat.Data[c*w.mat.Inc] = blas64.Dot(ar,
					blas64.Vector{Inc: 1, Data: a.Col(col, c)},
					b.mat,
				)
			}
		} else {
			row := make([]float64, ac)
			for r := 0; r < ar; r++ {
				w.mat.Data[r*w.mat.Inc] = blas64.Dot(ac,
					blas64.Vector{Inc: 1, Data: a.Row(row, r)},
					b.mat,
				)
			}
		}
	default:
		if trans {
			col := make([]float64, ar)
			for c := 0; c < ac; c++ {
				for i := range col {
					col[i] = a.At(i, c)
				}
				var f float64
				for i, e := range col {
					f += e * b.mat.Data[i*b.mat.Inc]
				}
				w.mat.Data[c*w.mat.Inc] = f
			}
		} else {
			row := make([]float64, ac)
			for r := 0; r < ar; r++ {
				for i := range row {
					row[i] = a.At(r, i)
				}
//.........这里部分代码省略.........

// svdCheck checks that the singular value decomposition correctly multiplies back
// to the original matrix.
func svdCheck(t *testing.T, thin bool, errStr string, m, n int, s, a, u []float64, ldu int, vt []float64, ldvt int, aCopy []float64, lda int) {
	sigma := blas64.General{
		Rows:   m,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, m*n),
	}
	for i := 0; i < min(m, n); i++ {
		sigma.Data[i*sigma.Stride+i] = s[i]
	}

	uMat := blas64.General{
		Rows:   m,
		Cols:   m,
		Stride: ldu,
		Data:   u,
	}
	vTMat := blas64.General{
		Rows:   n,
		Cols:   n,
		Stride: ldvt,
		Data:   vt,
	}
	if thin {
		sigma.Rows = min(m, n)
		sigma.Cols = min(m, n)
		uMat.Cols = min(m, n)
		vTMat.Rows = min(m, n)
	}

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

	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, uMat, sigma, 0, tmp)
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp, vTMat, 0, ans)

	if !floats.EqualApprox(ans.Data, aCopy, 1e-8) {
		t.Errorf("Decomposition mismatch. Trim = %v, %s", thin, errStr)
	}

	if !thin {
		// Check that U and V are orthogonal.
		for i := 0; i < uMat.Rows; i++ {
			for j := i + 1; j < uMat.Rows; j++ {
				dot := blas64.Dot(uMat.Cols,
					blas64.Vector{Inc: 1, Data: uMat.Data[i*uMat.Stride:]},
					blas64.Vector{Inc: 1, Data: uMat.Data[j*uMat.Stride:]},
				)
				if dot > 1e-8 {
					t.Errorf("U not orthogonal %s", errStr)
				}
			}
		}
		for i := 0; i < vTMat.Rows; i++ {
			for j := i + 1; j < vTMat.Rows; j++ {
				dot := blas64.Dot(vTMat.Cols,
					blas64.Vector{Inc: 1, Data: vTMat.Data[i*vTMat.Stride:]},
					blas64.Vector{Inc: 1, Data: vTMat.Data[j*vTMat.Stride:]},
				)
				if dot > 1e-8 {
					t.Errorf("V not orthogonal %s", errStr)
				}
			}
		}
	}
}

func Dgeqr2Test(t *testing.T, impl Dgeqr2er) {
	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},
		{12, 1, 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},
		{12, 1, 20},
	} {
		n := test.n
		m := test.m
		lda := test.lda
		if lda == 0 {
			lda = test.n
		}
		a := make([]float64, m*lda)
		for i := range a {
			a[i] = rand.Float64()
		}
		aCopy := make([]float64, len(a))
		k := min(m, n)
		tau := make([]float64, k)
		for i := range tau {
			tau[i] = rand.Float64()
		}
		work := make([]float64, n)
		for i := range work {
			work[i] = rand.Float64()
		}
		copy(aCopy, a)
		impl.Dgeqr2(m, n, a, lda, tau, work)

		// Test that the QR factorization has completed successfully. Compute
		// Q based on the vectors.
		q := constructQ("QR", m, n, a, lda, tau)

		// Check that q is orthonormal
		for i := 0; i < m; i++ {
			nrm := blas64.Nrm2(m, blas64.Vector{1, q.Data[i*m:]})
			if math.Abs(nrm-1) > 1e-14 {
				t.Errorf("Case %v, q not normal", c)
			}
			for j := 0; j < i; j++ {
				dot := blas64.Dot(m, blas64.Vector{1, q.Data[i*m:]}, blas64.Vector{1, q.Data[j*m:]})
				if math.Abs(dot) > 1e-14 {
					t.Errorf("Case %v, q not orthogonal", i)
				}
			}
		}
		// Check that A = Q * R
		r := blas64.General{
			Rows:   m,
			Cols:   n,
			Stride: n,
			Data:   make([]float64, m*n),
		}
		for i := 0; i < m; i++ {
			for j := i; j < n; j++ {
				r.Data[i*n+j] = a[i*lda+j]
			}
		}
		atmp := blas64.General{
			Rows:   m,
			Cols:   n,
			Stride: lda,
			Data:   make([]float64, m*lda),
		}
		copy(atmp.Data, a)
		blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, r, 0, atmp)
		if !floats.EqualApprox(atmp.Data, aCopy, 1e-14) {
			t.Errorf("Q*R != a")
		}
	}
}

//.........这里部分代码省略.........
			blas64.Trmm(blas.Left, aT, 1, amat, m.mat)
			return
		}
		if aU, ok := aU.(*Vector); ok {
			avec := aU.RawVector()
			if aTrans {
				// {1,ac} x {ac, bc}
				// Transpose B so that the vector is on the right.
				cvec := blas64.Vector{
					Inc:  1,
					Data: m.mat.Data,
				}
				bT := blas.Trans
				if bTrans {
					bT = blas.NoTrans
				}
				blas64.Gemv(bT, 1, bmat, avec, 0, cvec)
				return
			}
			// {ar,1} x {1,bc} which is not a vector result.
			// Instead, construct A as a General.
			amat := blas64.General{
				Rows:   ar,
				Cols:   1,
				Stride: avec.Inc,
				Data:   avec.Data,
			}
			blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
			return
		}
	}

	if aU, ok := aU.(Vectorer); ok {
		if bU, ok := bU.(Vectorer); ok {
			row := make([]float64, ac)
			col := make([]float64, br)
			if aTrans {
				if bTrans {
					for r := 0; r < ar; r++ {
						dataTmp := m.mat.Data[r*m.mat.Stride : r*m.mat.Stride+bc]
						for c := 0; c < bc; c++ {
							dataTmp[c] = blas64.Dot(ac,
								blas64.Vector{Inc: 1, Data: aU.Col(row, r)},
								blas64.Vector{Inc: 1, Data: bU.Row(col, c)},
							)
						}
					}
					return
				}
				// TODO(jonlawlor): determine if (b*a)' is more efficient
				for r := 0; r < ar; r++ {
					dataTmp := m.mat.Data[r*m.mat.Stride : r*m.mat.Stride+bc]
					for c := 0; c < bc; c++ {
						dataTmp[c] = blas64.Dot(ac,
							blas64.Vector{Inc: 1, Data: aU.Col(row, r)},
							blas64.Vector{Inc: 1, Data: bU.Col(col, c)},
						)
					}
				}
				return
			}
			if bTrans {
				for r := 0; r < ar; r++ {
					dataTmp := m.mat.Data[r*m.mat.Stride : r*m.mat.Stride+bc]
					for c := 0; c < bc; c++ {
						dataTmp[c] = blas64.Dot(ac,
							blas64.Vector{Inc: 1, Data: aU.Row(row, r)},
							blas64.Vector{Inc: 1, Data: bU.Row(col, c)},
						)
					}
				}
				return
			}
			for r := 0; r < ar; r++ {
				dataTmp := m.mat.Data[r*m.mat.Stride : r*m.mat.Stride+bc]
				for c := 0; c < bc; c++ {
					dataTmp[c] = blas64.Dot(ac,
						blas64.Vector{Inc: 1, Data: aU.Row(row, r)},
						blas64.Vector{Inc: 1, Data: bU.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)
			}
			m.mat.Data[r*m.mat.Stride+c] = v
		}
	}
}

//.........这里部分代码省略.........
				if uplo == blas.Upper {
					vi = blas64.Vector{
						Inc:  1,
						Data: make([]float64, n),
					}
					for j := 0; j < i; j++ {
						vi.Data[j] = a[j*lda+i+1]
					}
					vi.Data[i] = 1
				} else {
					vi = blas64.Vector{
						Inc:  1,
						Data: make([]float64, n),
					}
					vi.Data[i+1] = 1
					for j := i + 2; j < n; j++ {
						vi.Data[j] = a[j*lda+i]
					}
				}
				blas64.Ger(-tau[i], vi, vi, hMat)
				copy(qCopy.Data, qMat.Data)

				// Multiply q by the new h.
				if uplo == blas.Upper {
					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hMat, qCopy, 0, qMat)
				} else {
					blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, hMat, 0, qMat)
				}
			}
			// Check that Q is orthonormal
			othonormal := true
			for i := 0; i < n; i++ {
				for j := i; j < n; j++ {
					dot := blas64.Dot(n,
						blas64.Vector{Inc: 1, Data: qMat.Data[i*qMat.Stride:]},
						blas64.Vector{Inc: 1, Data: qMat.Data[j*qMat.Stride:]},
					)
					if i == j {
						if math.Abs(dot-1) > 1e-10 {
							othonormal = false
						}
					} else {
						if math.Abs(dot) > 1e-10 {
							othonormal = false
						}
					}
				}
			}
			if !othonormal {
				t.Errorf("Q not orthonormal")
			}

			// Compute Q^T * A * Q.
			aMat := blas64.General{
				Rows:   n,
				Cols:   n,
				Stride: n,
				Data:   make([]float64, len(a)),
			}

			for i := 0; i < n; i++ {
				for j := i; j < n; j++ {
					v := aCopy[i*lda+j]
					if uplo == blas.Lower {
						v = aCopy[j*lda+i]
					}

// MulTrans takes the matrix product of a and b, optionally transposing each,
// and placing the result in the receiver.
//
// See the MulTranser interface for more information.
func (m *Dense) MulTrans(a Matrix, aTrans bool, b Matrix, bTrans bool) {
	ar, ac := a.Dims()
	if aTrans {
		ar, ac = ac, ar
	}

	br, bc := b.Dims()
	if bTrans {
		br, bc = bc, br
	}

	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 := a.RawMatrix()
			if a == b && aTrans != bTrans {
				var op blas.Transpose
				if aTrans {
					op = blas.Trans
				} else {
					op = blas.NoTrans
				}
				blas64.Syrk(op, 1, amat, 0, blas64.Symmetric{N: w.Mat.Rows, Stride: w.Mat.Stride, Data: w.Mat.Data, Uplo: blas.Upper})

				// Fill lower matrix with result.
				// TODO(kortschak): Investigate whether using blas64.Copy improves the performance of this significantly.
				for i := 0; i < w.Mat.Rows; i++ {
					for j := i + 1; j < w.Mat.Cols; j++ {
						w.set(j, i, w.at(i, j))
					}
				}
			} else {
				var aOp, bOp blas.Transpose
				if aTrans {
					aOp = blas.Trans
				} else {
					aOp = blas.NoTrans
				}
				if bTrans {
					bOp = blas.Trans
				} else {
					bOp = blas.NoTrans
				}
				bmat := b.RawMatrix()
				blas64.Gemm(aOp, bOp, 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)
			if aTrans {
				if bTrans {
					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.Col(row, r)},
								blas64.Vector{Inc: 1, Data: b.Row(col, c)},
							)
						}
					}
					return
				}
				// TODO(jonlawlor): determine if (b*a)' is more efficient
				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.Col(row, r)},
							blas64.Vector{Inc: 1, Data: b.Col(col, c)},
						)
					}
				}
				return
			}
			if bTrans {
				for r := 0; r < ar; r++ {
//.........这里部分代码省略.........

// Mul multiplies two matrix and saves the result in m. Note that the
// arguments a or b should be either Matrix or *Dense.
// Therfore, if a or b is of type Dense, you'll need to pass them by address.
// For example: m.Mul(a, &b) when a is *Dense and b is Dense.
func (m *Dense) Mul(a, b Matrix) {
	ar, ac := a.Dims()
	br, bc := b.Dims()

	if ac != br {
		panic(ErrShape)
	}

	var w Dense
	if m != a && m != b {
		w = *m
	}
	if w.isZero() {
		w.mat = blas64.General{
			Rows:   ar,
			Cols:   bc,
			Stride: bc,
			Data:   use(w.mat.Data, ar*bc),
		}
	} else if ar != w.mat.Rows || bc != w.mat.Cols {
		panic(ErrShape)
	}

	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)
			*m = w
			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)},
					)
				}
			}
			*m = w
			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
		}
	}
	*m = w
}