text: // Copyright ©2016 The Gonum Authors. All rights reserved.
  // Use of this source code is governed by a BSD-style
  // license that can be found in the LICENSE file.
  
  package gonum
  
  import (
  	"math"
  
  	"gonum.org/v1/gonum/blas"
  	"gonum.org/v1/gonum/blas/blas64"
  	"gonum.org/v1/gonum/lapack"
  )
  
  // Dsyev computes all eigenvalues and, optionally, the eigenvectors of a real
  // symmetric matrix A.
  //
  // w contains the eigenvalues in ascending order upon return. w must have length
  // at least n, and Dsyev will panic otherwise.
  //
  // On entry, a contains the elements of the symmetric matrix A in the triangular
  // portion specified by uplo. If jobz == lapack.EVCompute, a contains the
  // orthonormal eigenvectors of A on exit, otherwise jobz must be lapack.EVNone
  // and on exit the specified triangular region is overwritten.
  //
  // work is temporary storage, and lwork specifies the usable memory length. At minimum,
  // lwork >= 3*n-1, and Dsyev will panic otherwise. The amount of blocking is
  // limited by the usable length. If lwork == -1, instead of computing Dsyev the
  // optimal work length is stored into work[0].
  func (impl Implementation) Dsyev(jobz lapack.EVJob, uplo blas.Uplo, n int, a []float64, lda int, w, work []float64, lwork int) (ok bool) {
  	switch {
  	case jobz != lapack.EVNone && jobz != lapack.EVCompute:
  		panic(badEVJob)
  	case uplo != blas.Upper && uplo != blas.Lower:
  		panic(badUplo)
  	case n < 0:
  		panic(nLT0)
  	case lda < max(1, n):
  		panic(badLdA)
  	case lwork < max(1, 3*n-1) && lwork != -1:
  		panic(badLWork)
  	case len(work) < max(1, lwork):
  		panic(shortWork)
  	}
  
  	// Quick return if possible.
  	if n == 0 {
  		return true
  	}
  
  	var opts string
  	if uplo == blas.Upper {
  		opts = "U"
  	} else {
  		opts = "L"
  	}
  	nb := impl.Ilaenv(1, "DSYTRD", opts, n, -1, -1, -1)
  	lworkopt := max(1, (nb+2)*n)
  	if lwork == -1 {
  		work[0] = float64(lworkopt)
  		return
  	}
  
  	switch {
  	case len(a) < (n-1)*lda+n:
  		panic(shortA)
  	case len(w) < n:
  		panic(shortW)
  	}
  
  	if n == 1 {
  		w[0] = a[0]
  		work[0] = 2
  		if jobz == lapack.EVCompute {
  			a[0] = 1
  		}
  		return true
  	}
  
  	safmin := dlamchS
  	eps := dlamchP
  	smlnum := safmin / eps
  	bignum := 1 / smlnum
  	rmin := math.Sqrt(smlnum)
  	rmax := math.Sqrt(bignum)
  
  	// Scale matrix to allowable range, if necessary.
  	anrm := impl.Dlansy(lapack.MaxAbs, uplo, n, a, lda, work)
  	scaled := false
  	var sigma float64
  	if anrm > 0 && anrm < rmin {
  		scaled = true
  		sigma = rmin / anrm
  	} else if anrm > rmax {
  		scaled = true
  		sigma = rmax / anrm
  	}
  	if scaled {
  		kind := lapack.LowerTri
  		if uplo == blas.Upper {
  			kind = lapack.UpperTri
  		}
  		impl.Dlascl(kind, 0, 0, 1, sigma, n, n, a, lda)
  	}
  	var inde int
  	indtau := inde + n
  	indwork := indtau + n
  	llwork := lwork - indwork
  	impl.Dsytrd(uplo, n, a, lda, w, work[inde:], work[indtau:], work[indwork:], llwork)
  
  	// For eigenvalues only, call Dsterf. For eigenvectors, first call Dorgtr
  	// to generate the orthogonal matrix, then call Dsteqr.
  	if jobz == lapack.EVNone {
  		ok = impl.Dsterf(n, w, work[inde:])
  	} else {
  		impl.Dorgtr(uplo, n, a, lda, work[indtau:], work[indwork:], llwork)
  		ok = impl.Dsteqr(lapack.EVComp(jobz), n, w, work[inde:], a, lda, work[indtau:])
  	}
  	if !ok {
  		return false
  	}
  
  	// If the matrix was scaled, then rescale eigenvalues appropriately.
  	if scaled {
  		bi := blas64.Implementation()
  		bi.Dscal(n, 1/sigma, w, 1)
  	}
  	work[0] = float64(lworkopt)
  	return true
  }
  