func (wm *writtenMemory) backwardWErase() {
	n := wm.N
	m := len(wm.TopVal) / n

	mgrad := make([]float64, n*m)
	mGradG := blas64.General{Rows: n, Cols: m, Stride: m, Data: mgrad}
	hEraseGrad := blas64.Vector{Inc: 1, Data: make([]float64, m)}
	for i, weights := range wm.Ws {
		erase := wm.erase[i]
		add := wm.add[i]
		eraseV := blas64.Vector{Inc: 1, Data: erase}
		addV := blas64.Vector{Inc: 1, Data: add}
		weightsVal := blas64.Vector{Inc: 1, Data: weights.TopVal}

		for j := range mgrad {
			mgrad[j] = 0
		}
		blas64.Ger(1, weightsVal, eraseV, mGradG)
		wm.div1MWE(mgrad)
		floats.Mul(mgrad, wm.TopGrad)

		weightsV := blas64.Vector{Inc: 1, Data: weights.TopGrad}
		blas64.Gemv(blas.NoTrans, -1, mGradG, eraseV, 1, weightsV)
		blas64.Gemv(blas.NoTrans, 1, blas64.General{Rows: n, Cols: m, Stride: m, Data: wm.TopGrad}, addV, 1, weightsV)

		hErase := wm.Heads[i].EraseGrad()
		for j := range hEraseGrad.Data {
			hEraseGrad.Data[j] = 0
		}
		blas64.Gemv(blas.Trans, -1, mGradG, weightsVal, 1, hEraseGrad)
		for j, e := range erase {
			hErase[j] += hEraseGrad.Data[j] * e * (1 - e)
		}
	}
}

// isLeftEigenvectorOf returns whether the vector yRe+i*yIm, where i is the
// imaginary unit, is the left eigenvector of A corresponding to the eigenvalue
// lambda.
//
// A left eigenvector corresponding to a complex eigenvalue λ is a complex
// non-zero vector y such that
//  y^H A = λ y^H,
// which is equivalent for real A to
//  A^T y = conj(λ) y,
func isLeftEigenvectorOf(a blas64.General, yRe, yIm []float64, lambda complex128, tol float64) bool {
	if a.Rows != a.Cols {
		panic("matrix not square")
	}

	if imag(lambda) != 0 && yIm == nil {
		// Complex eigenvalue of a real matrix cannot have a real
		// eigenvector.
		return false
	}

	n := a.Rows

	// Compute A^T real(y) and store the result into yReAns.
	yReAns := make([]float64, n)
	blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yRe}, 0, blas64.Vector{1, yReAns})

	if imag(lambda) == 0 && yIm == nil {
		// Real eigenvalue and eigenvector.

		// Compute λy and store the result into lambday.
		lambday := make([]float64, n)
		floats.AddScaled(lambday, real(lambda), yRe)

		if floats.Distance(yReAns, lambday, math.Inf(1)) > tol {
			return false
		}
		return true
	}

	// Complex eigenvector, and real or complex eigenvalue.

	// Compute A^T imag(y) and store the result into yImAns.
	yImAns := make([]float64, n)
	blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yIm}, 0, blas64.Vector{1, yImAns})

	// Compute conj(λ)y and store the result into lambday.
	lambda = cmplx.Conj(lambda)
	lambday := make([]complex128, n)
	for i := range lambday {
		lambday[i] = lambda * complex(yRe[i], yIm[i])
	}

	for i, v := range lambday {
		ay := complex(yReAns[i], yImAns[i])
		if cmplx.Abs(v-ay) > tol {
			return false
		}
	}
	return true
}

func (old *controller1) Forward(reads []*memRead, x []float64) Controller {
	c := controller1{
		weightsVal:  old.weightsVal,
		weightsGrad: old.weightsGrad,
		Reads:       reads,
		X:           x,
		H1Val:       make([]float64, old.h1Size+1),
		H1Grad:      make([]float64, old.h1Size+1),
		heads:       make([]*Head, len(reads)),
		outVal:      make([]float64, old.wyRows()),
		outGrad:     make([]float64, old.wyRows()),

		numHeads: old.numHeads,
		memoryM:  old.memoryM,
		memoryN:  old.memoryN,
		xSize:    old.xSize,
		h1Size:   old.h1Size,
		ySize:    old.ySize,
	}

	ud := make([]float64, c.wh1Cols())
	for i, read := range reads {
		copy(ud[i*c.memoryM:], read.TopVal)
	}
	copy(ud[c.numHeads*c.memoryM:], c.X)
	ud[c.numHeads*c.memoryM+c.xSize] = 1
	c.ReadsXVal = blas64.Vector{Inc: 1, Data: ud}

	h1 := blas64.Vector{Inc: 1, Data: c.H1Val[0:c.h1Size]}
	blas64.Gemv(blas.NoTrans, 1, c.wh1Val(), c.ReadsXVal, 1, h1)
	for i, h := range c.H1Val[0:c.h1Size] {
		c.H1Val[i] = Sigmoid(h)
	}

	c.H1Val[c.h1Size] = 1
	h1 = blas64.Vector{Inc: 1, Data: c.H1Val}
	outV := blas64.Vector{Inc: 1, Data: c.outVal}
	blas64.Gemv(blas.NoTrans, 1, c.wyVal(), h1, 1, outV)

	hul := headUnitsLen(c.memoryM)
	for i := range c.heads {
		head := NewHead(c.memoryM)
		c.heads[i] = head
		start := c.ySize + i*hul
		head.vals = c.outVal[start : start+hul]
		head.grads = c.outGrad[start : start+hul]
	}

	return &c
}

// isRightEigenvectorOf returns whether the vector xRe+i*xIm, where i is the
// imaginary unit, is the right eigenvector of A corresponding to the eigenvalue
// lambda.
//
// A right eigenvector corresponding to a complex eigenvalue λ is a complex
// non-zero vector x such that
//  A x = λ x.
func isRightEigenvectorOf(a blas64.General, xRe, xIm []float64, lambda complex128, tol float64) bool {
	if a.Rows != a.Cols {
		panic("matrix not square")
	}

	if imag(lambda) != 0 && xIm == nil {
		// Complex eigenvalue of a real matrix cannot have a real
		// eigenvector.
		return false
	}

	n := a.Rows

	// Compute A real(x) and store the result into xReAns.
	xReAns := make([]float64, n)
	blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xRe}, 0, blas64.Vector{1, xReAns})

	if imag(lambda) == 0 && xIm == nil {
		// Real eigenvalue and eigenvector.

		// Compute λx and store the result into lambdax.
		lambdax := make([]float64, n)
		floats.AddScaled(lambdax, real(lambda), xRe)

		if floats.Distance(xReAns, lambdax, math.Inf(1)) > tol {
			return false
		}
		return true
	}

	// Complex eigenvector, and real or complex eigenvalue.

	// Compute A imag(x) and store the result into xImAns.
	xImAns := make([]float64, n)
	blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xIm}, 0, blas64.Vector{1, xImAns})

	// Compute λx and store the result into lambdax.
	lambdax := make([]complex128, n)
	for i := range lambdax {
		lambdax[i] = lambda * complex(xRe[i], xIm[i])
	}

	for i, v := range lambdax {
		ax := complex(xReAns[i], xImAns[i])
		if cmplx.Abs(v-ax) > tol {
			return false
		}
	}
	return true
}

// replaces x with Q.x
func (f LQFactor) applyQTo(x *Dense, trans bool) {
	nh, nc := f.LQ.Dims()
	m, n := x.Dims()
	if m != nc {
		panic(ErrShape)
	}
	proj := make([]float64, n)

	if trans {
		for k := nh - 1; k >= 0; k-- {
			hh := f.LQ.RawRowView(k)[k:]

			sub := x.View(k, 0, m-k, n).(*Dense)

			blas64.Gemv(blas.Trans,
				1, sub.Mat, blas64.Vector{Inc: 1, Data: hh},
				0, blas64.Vector{Inc: 1, Data: proj},
			)
			for i := k; i < m; i++ {
				row := x.RawRowView(i)
				blas64.Axpy(n, -hh[i-k],
					blas64.Vector{Inc: 1, Data: proj},
					blas64.Vector{Inc: 1, Data: row},
				)
			}
		}
	} else {
		for k := 0; k < nh; k++ {
			hh := f.LQ.RawRowView(k)[k:]

			sub := x.View(k, 0, m-k, n).(*Dense)

			blas64.Gemv(blas.Trans,
				1, sub.Mat, blas64.Vector{Inc: 1, Data: hh},
				0, blas64.Vector{Inc: 1, Data: proj},
			)
			for i := k; i < m; i++ {
				row := x.RawRowView(i)
				blas64.Axpy(n, -hh[i-k],
					blas64.Vector{Inc: 1, Data: proj},
					blas64.Vector{Inc: 1, Data: row},
				)
			}
		}
	}
}

func (c *controller1) Backward() {
	out := blas64.Vector{Inc: 1, Data: c.outGrad}
	h1Val := blas64.Vector{Inc: 1, Data: c.H1Val}
	h1Grad := blas64.Vector{Inc: 1, Data: c.H1Grad}
	blas64.Gemv(blas.Trans, 1, c.wyVal(), out, 1, h1Grad)
	blas64.Ger(1, out, h1Val, c.wyGrad())

	h1Val = blas64.Vector{Inc: 1, Data: c.H1Val[0:c.h1Size]}
	h1Grad = blas64.Vector{Inc: 1, Data: c.H1Grad[0:c.h1Size]}
	for i, v := range h1Val.Data {
		h1Grad.Data[i] *= v * (1 - v)
	}

	u := blas64.Vector{Inc: 1, Data: make([]float64, c.wh1Cols())}
	blas64.Gemv(blas.Trans, 1, c.wh1Val(), h1Grad, 1, u)
	blas64.Ger(1, h1Grad, c.ReadsXVal, c.wh1Grad())

	for i, read := range c.Reads {
		copy(read.TopGrad, u.Data[i*c.memoryM:(i+1)*c.memoryM])
	}
}

func (r *memRead) Backward() {
	n := r.Memory.N
	m := len(r.Memory.TopVal) / n

	grad := blas64.Vector{Inc: 1, Data: r.TopGrad}
	memVal := blas64.General{Rows: n, Cols: m, Stride: m, Data: r.Memory.TopVal}
	weightsGrad := blas64.Vector{Inc: 1, Data: r.W.TopGrad}
	blas64.Gemv(blas.NoTrans, 1, memVal, grad, 1, weightsGrad)

	memGrad := blas64.General{Rows: n, Cols: m, Stride: m, Data: r.Memory.TopGrad}
	weights := blas64.Vector{Inc: 1, Data: r.W.TopVal}
	blas64.Ger(1, weights, grad, memGrad)
}

// BestClassifier returns the classifier in the static
// pool whose output correlates the most highly with
// the given weight vector, as measured by absolute
// cosine distance.
//
// The list argument is ignored, since a StaticPool
// always uses the set of samples it was given when
// it was initialized.
func (s *StaticPool) BestClassifier(list SampleList, weights linalg.Vector) Classifier {
	vec := blas64.Vector{
		Inc:  1,
		Data: weights,
	}
	output := blas64.Vector{
		Inc:  1,
		Data: make([]float64, len(s.classifiers)),
	}
	blas64.Gemv(blas.NoTrans, 1, s.outputMatrix, vec, 0, output)
	largest := blas64.Iamax(len(s.classifiers), output)
	return s.classifiers[largest]
}

func newMemRead(w *refocus, memory *writtenMemory) *memRead {
	m := len(memory.TopVal) / memory.N
	r := memRead{
		W:       w,
		Memory:  memory,
		TopVal:  make([]float64, m),
		TopGrad: make([]float64, m),
	}

	weights := blas64.Vector{Inc: 1, Data: w.TopVal}
	mem := blas64.General{Rows: memory.N, Cols: m, Stride: m, Data: memory.TopVal}
	top := blas64.Vector{Inc: 1, Data: r.TopVal}
	blas64.Gemv(blas.Trans, 1, mem, weights, 1, top)

	return &r
}

// 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)
				}
//.........这里部分代码省略.........

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

	aU, aTrans := untranspose(a)
	bU, bTrans := untranspose(b)
	m.reuseAs(ar, bc)
	var restore func()
	if m == aU {
		m, restore = m.isolatedWorkspace(aU)
		defer restore()
	} else if m == bU {
		m, restore = m.isolatedWorkspace(bU)
		defer restore()
	}
	aT := blas.NoTrans
	if aTrans {
		aT = blas.Trans
	}
	bT := blas.NoTrans
	if bTrans {
		bT = blas.Trans
	}

	// Some of the cases do not have a transpose option, so create
	// temporary memory.
	// C = A^T * B = (B^T * A)^T
	// C^T = B^T * A.
	if aU, ok := aU.(RawMatrixer); ok {
		amat := aU.RawMatrix()
		if bU, ok := bU.(RawMatrixer); ok {
			bmat := bU.RawMatrix()
			blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
			return
		}
		if bU, ok := bU.(RawSymmetricer); ok {
			bmat := bU.RawSymmetric()
			if aTrans {
				c := getWorkspace(ac, ar, false)
				blas64.Symm(blas.Left, 1, bmat, amat, 0, c.mat)
				strictCopy(m, c.T())
				putWorkspace(c)
				return
			}
			blas64.Symm(blas.Right, 1, bmat, amat, 0, m.mat)
			return
		}
		if bU, ok := bU.(RawTriangular); ok {
			// Trmm updates in place, so copy aU first.
			bmat := bU.RawTriangular()
			if aTrans {
				c := getWorkspace(ac, ar, false)
				var tmp Dense
				tmp.SetRawMatrix(aU.RawMatrix())
				c.Copy(&tmp)
				bT := blas.Trans
				if bTrans {
					bT = blas.NoTrans
				}
				blas64.Trmm(blas.Left, bT, 1, bmat, c.mat)
				strictCopy(m, c.T())
				putWorkspace(c)
				return
			}
			m.Copy(a)
			blas64.Trmm(blas.Right, bT, 1, bmat, m.mat)
			return
		}
		if bU, ok := bU.(*Vector); ok {
			bvec := bU.RawVector()
			if bTrans {
				// {ar,1} x {1,bc}, which is not a vector.
				// Instead, construct B as a General.
				bmat := blas64.General{
					Rows:   bc,
					Cols:   1,
					Stride: bvec.Inc,
					Data:   bvec.Data,
				}
				blas64.Gemm(aT, bT, 1, amat, bmat, 0, m.mat)
				return
			}
			cvec := blas64.Vector{
				Inc:  m.mat.Stride,
				Data: m.mat.Data,
			}
			blas64.Gemv(aT, 1, amat, bvec, 0, cvec)
			return
		}
	}
	if bU, ok := bU.(RawMatrixer); ok {
		bmat := bU.RawMatrix()
		if aU, ok := aU.(RawSymmetricer); ok {
//.........这里部分代码省略.........

// 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)
				}
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		{
			alpha: 0,
			n:     3,
		},
		{
			alpha: 1,
			n:     1,
		},
		{
			alpha: 1,
			n:     2,
			x:     []float64{4, 5, 6},
		},
	} {
		n := test.n
		incX := 1
		var x []float64
		if test.x == nil {
			x = make([]float64, n-1)
			for i := range x {
				x[i] = rand.Float64()
			}
		} else {
			x = make([]float64, n-1)
			copy(x, test.x)
		}
		xcopy := make([]float64, n-1)
		copy(xcopy, x)
		alpha := test.alpha
		beta, tau := impl.Dlarfg(n, alpha, x, incX)

		// Verify the returns and the values in v. Construct h and perform
		// the explicit multiplication.
		h := make([]float64, n*n)
		for i := 0; i < n; i++ {
			h[i*n+i] = 1
		}
		hmat := blas64.General{
			Rows:   n,
			Cols:   n,
			Stride: n,
			Data:   h,
		}
		v := make([]float64, n)
		copy(v[1:], x)
		v[0] = 1
		vVec := blas64.Vector{
			Inc:  1,
			Data: v,
		}
		blas64.Ger(-tau, vVec, vVec, hmat)
		eye := blas64.General{
			Rows:   n,
			Cols:   n,
			Stride: n,
			Data:   make([]float64, n*n),
		}
		blas64.Gemm(blas.Trans, blas.NoTrans, 1, hmat, hmat, 0, eye)
		iseye := true
		for i := 0; i < n; i++ {
			for j := 0; j < n; j++ {
				if i == j {
					if math.Abs(eye.Data[i*n+j]-1) > 1e-14 {
						iseye = false
					}
				} else {
					if math.Abs(eye.Data[i*n+j]) > 1e-14 {
						iseye = false
					}
				}
			}
		}
		if !iseye {
			t.Errorf("H^T * H is not I %V", eye)
		}

		xVec := blas64.Vector{
			Inc:  1,
			Data: make([]float64, n),
		}
		xVec.Data[0] = test.alpha
		copy(xVec.Data[1:], xcopy)

		ans := make([]float64, n)
		ansVec := blas64.Vector{
			Inc:  1,
			Data: ans,
		}
		blas64.Gemv(blas.NoTrans, 1, hmat, xVec, 0, ansVec)
		if math.Abs(ans[0]-beta) > 1e-14 {
			t.Errorf("Case %v, beta mismatch. Want %v, got %v", i, ans[0], beta)
		}
		for i := 1; i < n; i++ {
			if math.Abs(ans[i]) > 1e-14 {
				t.Errorf("Case %v, nonzero answer %v", i, ans)
				break
			}
		}
	}
}

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