// Row copies the elements in the jth column of the matrix into the slice dst.
// The length of the provided slice must equal the number of columns, unless the
// slice is nil in which case a new slice is first allocated.
func Row(dst []float64, i int, a Matrix) []float64 {
	r, c := a.Dims()
	if i < 0 || i >= r {
		panic(ErrColAccess)
	}
	if dst == nil {
		dst = make([]float64, c)
	} else {
		if len(dst) != c {
			panic(ErrColLength)
		}
	}
	aMat, aTrans := untranspose(a)
	if rm, ok := aMat.(RawMatrixer); ok {
		m := rm.RawMatrix()
		if aTrans {
			blas64.Copy(c,
				blas64.Vector{Inc: m.Stride, Data: m.Data[i:]},
				blas64.Vector{Inc: 1, Data: dst},
			)
			return dst
		}
		copy(dst, m.Data[i*m.Stride:i*m.Stride+m.Cols])
		return dst
	}
	for j := 0; j < c; j++ {
		dst[j] = a.At(i, j)
	}
	return dst
}

// CopyVec makes a copy of elements of a into the receiver. It is similar to the
// built-in copy; it copies as much as the overlap between the two vectors and
// returns the number of elements it copied.
func (v *Vector) CopyVec(a *Vector) int {
	n := min(v.Len(), a.Len())
	if v != a {
		blas64.Copy(n, a.mat, v.mat)
	}
	return n
}

// Col copies the elements in the jth column of the matrix into the slice dst.
// The length of the provided slice must equal the number of rows, unless the
// slice is nil in which case a new slice is first allocated.
func Col(dst []float64, j int, a Matrix) []float64 {
	r, c := a.Dims()
	if j < 0 || j >= c {
		panic(ErrColAccess)
	}
	if dst == nil {
		dst = make([]float64, r)
	} else {
		if len(dst) != r {
			panic(ErrRowLength)
		}
	}
	aMat, aTrans := untranspose(a)
	if rm, ok := aMat.(RawMatrixer); ok {
		m := rm.RawMatrix()
		if aTrans {
			copy(dst, m.Data[j*m.Stride:j*m.Stride+m.Cols])
			return dst
		}
		blas64.Copy(r,
			blas64.Vector{Inc: m.Stride, Data: m.Data[j:]},
			blas64.Vector{Inc: 1, Data: dst},
		)
		return dst
	}
	for i := 0; i < r; i++ {
		dst[i] = a.At(i, j)
	}
	return dst
}

// ScaleVec scales the vector a by alpha, placing the result in the receiver.
func (v *Vector) ScaleVec(alpha float64, a *Vector) {
	n := a.Len()
	if v != a {
		v.reuseAs(n)
		blas64.Copy(n, a.mat, v.mat)
	}
	blas64.Scal(n, alpha, v.mat)
}

// Clone makes a copy of a into the receiver, overwriting the previous value of
// the receiver. The clone operation does not make any restriction on shape and
// will not cause shadowing.
//
// See the Cloner interface for more information.
func (m *Dense) Clone(a Matrix) {
	r, c := a.Dims()
	mat := blas64.General{
		Rows:   r,
		Cols:   c,
		Stride: c,
	}
	m.capRows, m.capCols = r, c

	aU, trans := untranspose(a)
	switch aU := aU.(type) {
	case RawMatrixer:
		amat := aU.RawMatrix()
		mat.Data = make([]float64, r*c)
		if trans {
			for i := 0; i < r; i++ {
				blas64.Copy(c,
					blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
					blas64.Vector{Inc: 1, Data: mat.Data[i*c : (i+1)*c]})
			}
		} else {
			for i := 0; i < r; i++ {
				copy(mat.Data[i*c:(i+1)*c], amat.Data[i*amat.Stride:i*amat.Stride+c])
			}
		}
	case *Vector:
		amat := aU.mat
		mat.Data = make([]float64, aU.n)
		blas64.Copy(aU.n,
			blas64.Vector{Inc: amat.Inc, Data: amat.Data},
			blas64.Vector{Inc: 1, Data: mat.Data})
	default:
		mat.Data = make([]float64, r*c)
		w := *m
		w.mat = mat
		for i := 0; i < r; i++ {
			for j := 0; j < c; j++ {
				w.set(i, j, a.At(i, j))
			}
		}
		*m = w
		return
	}
	m.mat = mat
}

// CloneVec makes a copy of a into the receiver, overwriting the previous value
// of the receiver.
func (v *Vector) CloneVec(a *Vector) {
	if v == a {
		return
	}
	v.n = a.n
	v.mat = blas64.Vector{
		Inc:  1,
		Data: use(v.mat.Data, v.n),
	}
	blas64.Copy(v.n, a.mat, v.mat)
}

// Clone makes a copy of a into the receiver, overwriting the previous value of
// the receiver. The clone operation does not make any restriction on shape.
//
// See the Cloner interface for more information.
func (m *Dense) Clone(a Matrix) {
	r, c := a.Dims()
	mat := blas64.General{
		Rows:   r,
		Cols:   c,
		Stride: c,
	}
	m.capRows, m.capCols = r, c

	aU, trans := untranspose(a)
	switch aU := aU.(type) {
	case RawMatrixer:
		amat := aU.RawMatrix()
		// TODO(kortschak): Consider being more precise with determining whether a and m are aliases.
		// The current approach is that all RawMatrixers are considered potential aliases.
		// Note that below we assume that non-RawMatrixers are not aliases; this is not necessarily
		// true, but cases where it is not are not sensible. We should probably fix or document
		// this though.
		mat.Data = make([]float64, r*c)
		if trans {
			for i := 0; i < r; i++ {
				blas64.Copy(c,
					blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
					blas64.Vector{Inc: 1, Data: mat.Data[i*c : (i+1)*c]})
			}
		} else {
			for i := 0; i < r; i++ {
				copy(mat.Data[i*c:(i+1)*c], amat.Data[i*amat.Stride:i*amat.Stride+c])
			}
		}
	case Vectorer:
		mat.Data = use(m.mat.Data, r*c)
		if trans {
			for i := 0; i < r; i++ {
				aU.Col(mat.Data[i*c:(i+1)*c], i)
			}
		} else {
			for i := 0; i < r; i++ {
				aU.Row(mat.Data[i*c:(i+1)*c], i)
			}
		}
	default:
		mat.Data = use(m.mat.Data, r*c)
		m.mat = mat
		for i := 0; i < r; i++ {
			for j := 0; j < c; j++ {
				m.set(i, j, a.At(i, j))
			}
		}
		return
	}
	m.mat = mat
}

// SetCol sets the elements of the matrix in the specified column to the values
// of src.
//
// See the VectorSetter interface for more information.
func (m *Dense) SetCol(j int, src []float64) int {
	if j >= m.mat.Cols || j < 0 {
		panic(ErrColAccess)
	}

	blas64.Copy(min(len(src), m.mat.Rows),
		blas64.Vector{Inc: 1, Data: src},
		blas64.Vector{Inc: m.mat.Stride, Data: m.mat.Data[j:]},
	)

	return min(len(src), m.mat.Rows)
}

// SetCol sets the values in the specified column of the matrix to the values
// in src. len(src) must equal the number of rows in the receiver.
func (m *Dense) SetCol(j int, src []float64) {
	if j >= m.mat.Cols || j < 0 {
		panic(matrix.ErrColAccess)
	}
	if len(src) != m.mat.Rows {
		panic(matrix.ErrColLength)
	}

	blas64.Copy(m.mat.Rows,
		blas64.Vector{Inc: 1, Data: src},
		blas64.Vector{Inc: m.mat.Stride, Data: m.mat.Data[j:]},
	)
}

// AddScaledVec adds the vectors a and alpha*b, placing the result in the receiver.
func (v *Vector) AddScaledVec(a *Vector, alpha float64, b *Vector) {
	if alpha == 1 {
		v.AddVec(a, b)
		return
	}
	if alpha == -1 {
		v.SubVec(a, b)
		return
	}

	ar := a.Len()
	br := b.Len()

	if ar != br {
		panic(matrix.ErrShape)
	}

	v.reuseAs(ar)

	if alpha == 0 {
		v.CopyVec(a)
		return
	}

	switch {
	case v == a && v == b: // v <- v + alpha * v = (alpha + 1) * v
		blas64.Scal(ar, alpha+1, v.mat)
	case v == a && v != b: // v <- v + alpha * b
		blas64.Axpy(ar, alpha, b.mat, v.mat)
	case v != a && v == b: // v <- a + alpha * v
		if v.mat.Inc == 1 && a.mat.Inc == 1 {
			// Fast path for a common case.
			v := v.mat.Data
			for i, a := range a.mat.Data {
				v[i] *= alpha
				v[i] += a
			}
			return
		}
		blas64.Scal(ar, alpha, v.mat)
		blas64.Axpy(ar, 1, a.mat, v.mat)
	default: // v <- a + alpha * b
		if v.mat.Inc == 1 && a.mat.Inc == 1 && b.mat.Inc == 1 {
			// Fast path for a common case.
			asm.DaxpyUnitary(alpha, b.mat.Data, a.mat.Data, v.mat.Data)
			return
		}
		blas64.Copy(ar, a.mat, v.mat)
		blas64.Axpy(ar, alpha, b.mat, v.mat)
	}
}

// Col copies the elements in the jth column of the matrix into the slice dst.
// If the provided slice is nil, a new slice is first allocated.
//
// See the Vectorer interface for more information.
func (m *Dense) Col(dst []float64, j int) []float64 {
	if j >= m.mat.Cols || j < 0 {
		panic(ErrColAccess)
	}

	if dst == nil {
		dst = make([]float64, m.mat.Rows)
	}
	dst = dst[:min(len(dst), m.mat.Rows)]
	blas64.Copy(len(dst),
		blas64.Vector{Inc: m.mat.Stride, Data: m.mat.Data[j:]},
		blas64.Vector{Inc: 1, Data: dst},
	)

	return dst
}

// Copy makes a copy of elements of a into the receiver. It is similar to the
// built-in copy; it copies as much as the overlap between the two matrices and
// returns the number of rows and columns it copied.
//
// See the Copier interface for more information.
func (m *Dense) Copy(a Matrix) (r, c int) {
	r, c = a.Dims()
	if a == m {
		return r, c
	}
	r = min(r, m.mat.Rows)
	c = min(c, m.mat.Cols)
	if r == 0 || c == 0 {
		return 0, 0
	}

	aU, trans := untranspose(a)
	switch aU := aU.(type) {
	case RawMatrixer:
		amat := aU.RawMatrix()
		if trans {
			for i := 0; i < r; i++ {
				blas64.Copy(c,
					blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
					blas64.Vector{Inc: 1, Data: m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]})
			}
		} else {
			for i := 0; i < r; i++ {
				copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
			}
		}
	case Vectorer:
		if trans {
			for i := 0; i < r; i++ {
				aU.Col(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], i)
			}
		} else {
			for i := 0; i < r; i++ {
				aU.Row(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], i)
			}
		}
	default:
		for i := 0; i < r; i++ {
			for j := 0; j < c; j++ {
				m.set(i, j, a.At(i, j))
			}
		}
	}

	return r, c
}

// Copy makes a copy of elements of a into the receiver. It is similar to the
// built-in copy; it copies as much as the overlap between the two matrices and
// returns the number of rows and columns it copied. If a aliases the receiver
// and is a transposed Dense or Vector, with a non-unitary increment, Copy will
// panic.
//
// See the Copier interface for more information.
func (m *Dense) Copy(a Matrix) (r, c int) {
	r, c = a.Dims()
	if a == m {
		return r, c
	}
	r = min(r, m.mat.Rows)
	c = min(c, m.mat.Cols)
	if r == 0 || c == 0 {
		return 0, 0
	}

	aU, trans := untranspose(a)
	switch aU := aU.(type) {
	case RawMatrixer:
		amat := aU.RawMatrix()
		if trans {
			if amat.Stride != 1 {
				m.checkOverlap(amat)
			}
			for i := 0; i < r; i++ {
				blas64.Copy(c,
					blas64.Vector{Inc: amat.Stride, Data: amat.Data[i : i+(c-1)*amat.Stride+1]},
					blas64.Vector{Inc: 1, Data: m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]})
			}
		} else {
			switch o := offset(m.mat.Data, amat.Data); {
			case o < 0:
				for i := r - 1; i >= 0; i-- {
					copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
				}
			case o > 0:
				for i := 0; i < r; i++ {
					copy(m.mat.Data[i*m.mat.Stride:i*m.mat.Stride+c], amat.Data[i*amat.Stride:i*amat.Stride+c])
				}
			default:
				// Nothing to do.
			}
		}
	case *Vector:
		var n, stride int
		amat := aU.mat
		if trans {
			if amat.Inc != 1 {
				m.checkOverlap(aU.asGeneral())
			}
			n = c
			stride = 1
		} else {
			n = r
			stride = m.mat.Stride
		}
		if amat.Inc == 1 && stride == 1 {
			copy(m.mat.Data, amat.Data[:n])
			break
		}
		switch o := offset(m.mat.Data, amat.Data); {
		case o < 0:
			blas64.Copy(n,
				blas64.Vector{Inc: -amat.Inc, Data: amat.Data},
				blas64.Vector{Inc: -stride, Data: m.mat.Data})
		case o > 0:
			blas64.Copy(n,
				blas64.Vector{Inc: amat.Inc, Data: amat.Data},
				blas64.Vector{Inc: stride, Data: m.mat.Data})
		default:
			// Nothing to do.
		}
	default:
		for i := 0; i < r; i++ {
			for j := 0; j < c; j++ {
				m.set(i, j, a.At(i, j))
			}
		}
	}

	return r, c
}

// CopyVec makes a copy of elements of a into the receiver. It is similar to the
// built-in copy; it copies as much as the overlap between the two matrices and
// returns the number of rows and columns it copied.
func (v *Vector) CopyVec(a *Vector) (n int) {
	n = min(v.Len(), a.Len())
	blas64.Copy(n, a.mat, v.mat)
	return n
}

// SymRankOne performs a rank-1 update of the original matrix A and refactorizes
// its Cholesky factorization, storing the result into the reciever. That is, if
// in the original Cholesky factorization
//  U^T * U = A,
// in the updated factorization
//  U'^T * U' = A + alpha * x * x^T = A'.
//
// Note that when alpha is negative, the updating problem may be ill-conditioned
// and the results may be inaccurate, or the updated matrix A' may not be
// positive definite and not have a Cholesky factorization. SymRankOne returns
// whether the updated matrix A' is positive definite.
//
// SymRankOne updates a Cholesky factorization in O(n²) time. The Cholesky
// factorization computation from scratch is O(n³).
func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x *Vector) (ok bool) {
	if !orig.valid() {
		panic(badCholesky)
	}
	n := orig.Size()
	if x.Len() != n {
		panic(matrix.ErrShape)
	}
	if orig != c {
		if c.isZero() {
			c.chol = NewTriDense(n, matrix.Upper, nil)
		} else if c.chol.mat.N != n {
			panic(matrix.ErrShape)
		}
		c.chol.Copy(orig.chol)
	}

	if alpha == 0 {
		return true
	}

	// Algorithms for updating and downdating the Cholesky factorization are
	// described, for example, in
	// - J. J. Dongarra, J. R. Bunch, C. B. Moler, G. W. Stewart: LINPACK
	//   Users' Guide. SIAM (1979), pages 10.10--10.14
	// or
	// - P. E. Gill, G. H. Golub, W. Murray, and M. A. Saunders: Methods for
	//   modifying matrix factorizations. Mathematics of Computation 28(126)
	//   (1974), Method C3 on page 521
	//
	// The implementation is based on LINPACK code
	// http://www.netlib.org/linpack/dchud.f
	// http://www.netlib.org/linpack/dchdd.f
	// and
	// https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2646
	//
	// According to http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg00301.html
	// LINPACK is released under BSD license.
	//
	// See also:
	// - M. A. Saunders: Large-scale Linear Programming Using the Cholesky
	//   Factorization. Technical Report Stanford University (1972)
	//   http://i.stanford.edu/pub/cstr/reports/cs/tr/72/252/CS-TR-72-252.pdf
	// - Matthias Seeger: Low rank updates for the Cholesky decomposition.
	//   EPFL Technical Report 161468 (2004)
	//   http://infoscience.epfl.ch/record/161468

	work := make([]float64, n)
	blas64.Copy(n, x.RawVector(), blas64.Vector{1, work})

	if alpha > 0 {
		// Compute rank-1 update.
		if alpha != 1 {
			blas64.Scal(n, math.Sqrt(alpha), blas64.Vector{1, work})
		}
		umat := c.chol.mat
		stride := umat.Stride
		for i := 0; i < n; i++ {
			// Compute parameters of the Givens matrix that zeroes
			// the i-th element of x.
			c, s, r, _ := blas64.Rotg(umat.Data[i*stride+i], work[i])
			if r < 0 {
				// Multiply by -1 to have positive diagonal
				// elemnts.
				r *= -1
				c *= -1
				s *= -1
			}
			umat.Data[i*stride+i] = r
			if i < n-1 {
				// Multiply the extended factorization matrix by
				// the Givens matrix from the left. Only
				// the i-th row and x are modified.
				blas64.Rot(n-i-1,
					blas64.Vector{1, umat.Data[i*stride+i+1 : i*stride+n]},
					blas64.Vector{1, work[i+1 : n]},
					c, s)
			}
		}
		c.updateCond(-1)
		return true
	}

	// Compute rank-1 downdate.
	alpha = math.Sqrt(-alpha)
	if alpha != 1 {
//.........这里部分代码省略.........