vendor/gonum.org/v1/gonum/lapack/gonum/dsyev.go

   1 // Copyright ©2016 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package gonum
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12         "gonum.org/v1/gonum/lapack"
  13 )
  14
  15 // Dsyev computes all eigenvalues and, optionally, the eigenvectors of a real
  16 // symmetric matrix A.
  17 //
  18 // w contains the eigenvalues in ascending order upon return. w must have length
  19 // at least n, and Dsyev will panic otherwise.
  20 //
  21 // On entry, a contains the elements of the symmetric matrix A in the triangular
  22 // portion specified by uplo. If jobz == lapack.ComputeEV a contains the
  23 // orthonormal eigenvectors of A on exit, otherwise on exit the specified
  24 // triangular region is overwritten.
  25 //
  26 // work is temporary storage, and lwork specifies the usable memory length. At minimum,
  27 // lwork >= 3*n-1, and Dsyev will panic otherwise. The amount of blocking is
  28 // limited by the usable length. If lwork == -1, instead of computing Dsyev the
  29 // optimal work length is stored into work[0].
  30 func (impl Implementation) Dsyev(jobz lapack.EVJob, uplo blas.Uplo, n int, a []float64, lda int, w, work []float64, lwork int) (ok bool) {
  31         checkMatrix(n, n, a, lda)
  32         upper := uplo == blas.Upper
  33         wantz := jobz == lapack.ComputeEV
  34         var opts string
  35         if upper {
  36                 opts = "U"
  37         } else {
  38                 opts = "L"
  39         }
  40         nb := impl.Ilaenv(1, "DSYTRD", opts, n, -1, -1, -1)
  41         lworkopt := max(1, (nb+2)*n)
  42         work[0] = float64(lworkopt)
  43         if lwork == -1 {
  44                 return
  45         }
  46         if len(work) < lwork {
  47                 panic(badWork)
  48         }
  49         if lwork < 3*n-1 {
  50                 panic(badWork)
  51         }
  52         if n == 0 {
  53                 return true
  54         }
  55         if n == 1 {
  56                 w[0] = a[0]
  57                 work[0] = 2
  58                 if wantz {
  59                         a[0] = 1
  60                 }
  61                 return true
  62         }
  63         safmin := dlamchS
  64         eps := dlamchP
  65         smlnum := safmin / eps
  66         bignum := 1 / smlnum
  67         rmin := math.Sqrt(smlnum)
  68         rmax := math.Sqrt(bignum)
  69
  70         // Scale matrix to allowable range, if necessary.
  71         anrm := impl.Dlansy(lapack.MaxAbs, uplo, n, a, lda, work)
  72         scaled := false
  73         var sigma float64
  74         if anrm > 0 && anrm < rmin {
  75                 scaled = true
  76                 sigma = rmin / anrm
  77         } else if anrm > rmax {
  78                 scaled = true
  79                 sigma = rmax / anrm
  80         }
  81         if scaled {
  82                 kind := lapack.LowerTri
  83                 if upper {
  84                         kind = lapack.UpperTri
  85                 }
  86                 impl.Dlascl(kind, 0, 0, 1, sigma, n, n, a, lda)
  87         }
  88         var inde int
  89         indtau := inde + n
  90         indwork := indtau + n
  91         llwork := lwork - indwork
  92         impl.Dsytrd(uplo, n, a, lda, w, work[inde:], work[indtau:], work[indwork:], llwork)
  93
  94         // For eigenvalues only, call Dsterf. For eigenvectors, first call Dorgtr
  95         // to generate the orthogonal matrix, then call Dsteqr.
  96         if !wantz {
  97                 ok = impl.Dsterf(n, w, work[inde:])
  98         } else {
  99                 impl.Dorgtr(uplo, n, a, lda, work[indtau:], work[indwork:], llwork)
 100                 ok = impl.Dsteqr(lapack.EVComp(jobz), n, w, work[inde:], a, lda, work[indtau:])
 101         }
 102         if !ok {
 103                 return false
 104         }
 105
 106         // If the matrix was scaled, then rescale eigenvalues appropriately.
 107         if scaled {
 108                 bi := blas64.Implementation()
 109                 bi.Dscal(n, 1/sigma, w, 1)
 110         }
 111         work[0] = float64(lworkopt)
 112         return true
 113 }