func (s *SymDense) PowPSD(a Symmetric, pow float64) error {
	dim := a.Symmetric()
	s.reuseAs(dim)
	var eigen EigenSym
	ok := eigen.Factorize(a, true)
	if !ok {
		return ErrFailedEigen
	}
	values := eigen.Values(nil)
	for i, v := range values {
		if v <= 0 {
			return ErrNotPSD
		}
		values[i] = math.Pow(v, pow)
	}
	var u Dense
	u.EigenvectorsSym(&eigen)
	s.SymOuterK(values[0], u.ColView(0))
	var v VecDense
	for i := 1; i < dim; i++ {
		v.ColViewOf(&u, i)
		s.SymRankOne(s, values[i], &v)
	}
	return nil
}
PowPSD computes a^pow where a is a positive symmetric definite matrix.

PowPSD returns an error if the matrix is not  not positive symmetric definite
or the Eigendecomposition is not successful.
package mat
import (
	"math"
	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/blas/blas64"
)
var (
	symDense *SymDense
	_ Matrix           = symDense
	_ Symmetric        = symDense
	_ RawSymmetricer   = symDense
	_ MutableSymmetric = symDense
)
const (
	badSymTriangle = "mat: blas64.Symmetric not upper"
	badSymCap      = "mat: bad capacity for SymDense"
)
type SymDense struct {
	mat blas64.Symmetric
	cap int
}
type Symmetric interface {
	Matrix
	Symmetric() int
}
type RawSymmetricer interface {
	RawSymmetric() blas64.Symmetric
}
type MutableSymmetric interface {
	Symmetric
	SetSym(i, j int, v float64)
}
func NewSymDense(n int, data []float64) *SymDense {
	if n < 0 {
		panic("mat: negative dimension")
	}
	if data != nil && n*n != len(data) {
		panic(ErrShape)
	}
	if data == nil {
		data = make([]float64, n*n)
	}
	return &SymDense{
		mat: blas64.Symmetric{
			N:      n,
			Stride: n,
			Data:   data,
			Uplo:   blas.Upper,
		},
		cap: n,
	}
}
func (s *SymDense) Dims() (r, c int) {
	return s.mat.N, s.mat.N
}
func (s *SymDense) Caps() (r, c int) {
	return s.cap, s.cap
}
func (s *SymDense) T() Matrix {
	return s
}
func (s *SymDense) Symmetric() int {
	return s.mat.N
}
func (s *SymDense) RawSymmetric() blas64.Symmetric {
	return s.mat
}
func (s *SymDense) SetRawSymmetric(b blas64.Symmetric) {
	if b.Uplo != blas.Upper {
		panic(badSymTriangle)
	}
	s.mat = b
}
func (s *SymDense) Reset() {
	s.mat.N, s.mat.Stride = 0, 0
	s.mat.Data = s.mat.Data[:0]
}
func (s *SymDense) IsZero() bool {
	return s.mat.N == 0
}
func (s *SymDense) reuseAs(n int) {
	if n == 0 {
		panic(ErrZeroLength)
	}
	if s.mat.N > s.cap {
		panic(badSymCap)
	}
	if s.IsZero() {
		s.mat = blas64.Symmetric{
			N:      n,
			Stride: n,
			Data:   use(s.mat.Data, n*n),
			Uplo:   blas.Upper,
		}
		s.cap = n
		return
	}
	if s.mat.Uplo != blas.Upper {
		panic(badSymTriangle)
	}
	if s.mat.N != n {
		panic(ErrShape)
	}
}
func (s *SymDense) isolatedWorkspace(a Symmetric) (w *SymDense, restore func()) {
	n := a.Symmetric()
	if n == 0 {
		panic(ErrZeroLength)
	}
	w = getWorkspaceSym(n, false)
	return w, func() {
		s.CopySym(w)
		putWorkspaceSym(w)
	}
}
func (s *SymDense) AddSym(a, b Symmetric) {
	n := a.Symmetric()
	if n != b.Symmetric() {
		panic(ErrShape)
	}
	s.reuseAs(n)
	if a, ok := a.(RawSymmetricer); ok {
		if b, ok := b.(RawSymmetricer); ok {
			amat, bmat := a.RawSymmetric(), b.RawSymmetric()
			if s != a {
				s.checkOverlap(generalFromSymmetric(amat))
			}
			if s != b {
				s.checkOverlap(generalFromSymmetric(bmat))
			}
			for i := 0; i < n; i++ {
				btmp := bmat.Data[i*bmat.Stride+i : i*bmat.Stride+n]
				stmp := s.mat.Data[i*s.mat.Stride+i : i*s.mat.Stride+n]
				for j, v := range amat.Data[i*amat.Stride+i : i*amat.Stride+n] {
					stmp[j] = v + btmp[j]
				}
			}
			return
		}
	}
	s.checkOverlapMatrix(a)
	s.checkOverlapMatrix(b)
	for i := 0; i < n; i++ {
		stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
		for j := i; j < n; j++ {
			stmp[j] = a.At(i, j) + b.At(i, j)
		}
	}
}
func (s *SymDense) CopySym(a Symmetric) int {
	n := a.Symmetric()
	n = min(n, s.mat.N)
	if n == 0 {
		return 0
	}
	switch a := a.(type) {
	case RawSymmetricer:
		amat := a.RawSymmetric()
		if amat.Uplo != blas.Upper {
			panic(badSymTriangle)
		}
		for i := 0; i < n; i++ {
			copy(s.mat.Data[i*s.mat.Stride+i:i*s.mat.Stride+n], amat.Data[i*amat.Stride+i:i*amat.Stride+n])
		}
	default:
		for i := 0; i < n; i++ {
			stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
			for j := i; j < n; j++ {
				stmp[j] = a.At(i, j)
			}
		}
	}
	return n
}
func (s *SymDense) SymRankOne(a Symmetric, alpha float64, x Vector) {
	n, c := x.Dims()
	if a.Symmetric() != n || c != 1 {
		panic(ErrShape)
	}
	s.reuseAs(n)
	if s != a {
		if rs, ok := a.(RawSymmetricer); ok {
			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
		}
		s.CopySym(a)
	}
	xU, _ := untranspose(x)
	if rv, ok := xU.(RawVectorer); ok {
		xmat := rv.RawVector()
		s.checkOverlap((&VecDense{mat: xmat, n: n}).asGeneral())
		blas64.Syr(alpha, xmat, s.mat)
		return
	}
	for i := 0; i < n; i++ {
		for j := i; j < n; j++ {
			s.set(i, j, s.at(i, j)+alpha*x.AtVec(i)*x.AtVec(j))
		}
	}
}
func (s *SymDense) SymRankK(a Symmetric, alpha float64, x Matrix) {
	n := a.Symmetric()
	r, _ := x.Dims()
	if r != n {
		panic(ErrShape)
	}
	xMat, aTrans := untranspose(x)
	var g blas64.General
	if rm, ok := xMat.(RawMatrixer); ok {
		g = rm.RawMatrix()
	} else {
		g = DenseCopyOf(x).mat
		aTrans = false
	}
	if a != s {
		if rs, ok := a.(RawSymmetricer); ok {
			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
		}
		s.reuseAs(n)
		s.CopySym(a)
	}
	t := blas.NoTrans
	if aTrans {
		t = blas.Trans
	}
	blas64.Syrk(t, alpha, g, 1, s.mat)
}
func (s *SymDense) SymOuterK(alpha float64, x Matrix) {
	n, _ := x.Dims()
	switch {
	case s.IsZero():
		s.mat = blas64.Symmetric{
			N:      n,
			Stride: n,
			Data:   useZeroed(s.mat.Data, n*n),
			Uplo:   blas.Upper,
		}
		s.cap = n
		s.SymRankK(s, alpha, x)
	case s.mat.Uplo != blas.Upper:
		panic(badSymTriangle)
	case s.mat.N == n:
		if s == x {
			w := getWorkspaceSym(n, true)
			w.SymRankK(w, alpha, x)
			s.CopySym(w)
			putWorkspaceSym(w)
		} else {
			switch r := x.(type) {
			case RawMatrixer:
				s.checkOverlap(r.RawMatrix())
			case RawSymmetricer:
				s.checkOverlap(generalFromSymmetric(r.RawSymmetric()))
			case RawTriangular:
				s.checkOverlap(generalFromTriangular(r.RawTriangular()))
			}
			for i := 0; i < n; i++ {
				ri := i * s.mat.Stride
				zero(s.mat.Data[ri+i : ri+n])
			}
			s.SymRankK(s, alpha, x)
		}
	default:
		panic(ErrShape)
	}
}
func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
	n := s.mat.N
	xr, xc := x.Dims()
	if xr != n || xc != 1 {
		panic(ErrShape)
	}
	yr, yc := y.Dims()
	if yr != n || yc != 1 {
		panic(ErrShape)
	}
	if s != a {
		if rs, ok := a.(RawSymmetricer); ok {
			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
		}
	}
	var xmat, ymat blas64.Vector
	fast := true
	xU, _ := untranspose(x)
	if rv, ok := xU.(RawVectorer); ok {
		xmat = rv.RawVector()
		s.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
	} else {
		fast = false
	}
	yU, _ := untranspose(y)
	if rv, ok := yU.(RawVectorer); ok {
		ymat = rv.RawVector()
		s.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
	} else {
		fast = false
	}
	if s != a {
		if rs, ok := a.(RawSymmetricer); ok {
			s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
		}
		s.reuseAs(n)
		s.CopySym(a)
	}
	if fast {
		if s != a {
			s.reuseAs(n)
			s.CopySym(a)
		}
		blas64.Syr2(alpha, xmat, ymat, s.mat)
		return
	}
	for i := 0; i < n; i++ {
		s.reuseAs(n)
		for j := i; j < n; j++ {
			s.set(i, j, a.At(i, j)+alpha*(x.AtVec(i)*y.AtVec(j)+y.AtVec(i)*x.AtVec(j)))
		}
	}
}
func (s *SymDense) ScaleSym(f float64, a Symmetric) {
	n := a.Symmetric()
	s.reuseAs(n)
	if a, ok := a.(RawSymmetricer); ok {
		amat := a.RawSymmetric()
		if s != a {
			s.checkOverlap(generalFromSymmetric(amat))
		}
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				s.mat.Data[i*s.mat.Stride+j] = f * amat.Data[i*amat.Stride+j]
			}
		}
		return
	}
	for i := 0; i < n; i++ {
		for j := i; j < n; j++ {
			s.mat.Data[i*s.mat.Stride+j] = f * a.At(i, j)
		}
	}
}
func (s *SymDense) SubsetSym(a Symmetric, set []int) {
	n := len(set)
	na := a.Symmetric()
	s.reuseAs(n)
	var restore func()
	if a == s {
		s, restore = s.isolatedWorkspace(a)
		defer restore()
	}
	if a, ok := a.(RawSymmetricer); ok {
		raw := a.RawSymmetric()
		if s != a {
			s.checkOverlap(generalFromSymmetric(raw))
		}
		for i := 0; i < n; i++ {
			ssub := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
			r := set[i]
			rsub := raw.Data[r*raw.Stride : r*raw.Stride+na]
			for j := i; j < n; j++ {
				c := set[j]
				if r <= c {
					ssub[j] = rsub[c]
				} else {
					ssub[j] = raw.Data[c*raw.Stride+r]
				}
			}
		}
		return
	}
	for i := 0; i < n; i++ {
		for j := i; j < n; j++ {
			s.mat.Data[i*s.mat.Stride+j] = a.At(set[i], set[j])
		}
	}
}
func (s *SymDense) SliceSquare(i, k int) Matrix {
	sz := s.cap
	if i < 0 || sz < i || k < i || sz < k {
		panic(ErrIndexOutOfRange)
	}
	v := *s
	v.mat.Data = s.mat.Data[i*s.mat.Stride+i : (k-1)*s.mat.Stride+k]
	v.mat.N = k - i
	v.cap = s.cap - i
	return &v
}
func (s *SymDense) GrowSquare(n int) Matrix {
	if n < 0 {
		panic(ErrIndexOutOfRange)
	}
	if n == 0 {
		return s
	}
	var v SymDense
	n += s.mat.N
	if n > s.cap {
		v.mat = blas64.Symmetric{
			N:      n,
			Stride: n,
			Uplo:   blas.Upper,
			Data:   make([]float64, n*n),
		}
		v.cap = n
		var tmp SymDense
		tmp.mat = blas64.Symmetric{
			N:      s.cap,
			Stride: s.mat.Stride,
			Data:   s.mat.Data,
			Uplo:   s.mat.Uplo,
		}
		tmp.cap = s.cap
		v.CopySym(&tmp)
		return &v
	}
	v.mat = blas64.Symmetric{
		N:      n,
		Stride: s.mat.Stride,
		Uplo:   blas.Upper,
		Data:   s.mat.Data[:(n-1)*s.mat.Stride+n],
	}
	v.cap = s.cap
	return &v
}