github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dsyev.go (about)

     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  	"github.com/jingcheng-WU/gonum/blas"
    11  	"github.com/jingcheng-WU/gonum/blas/blas64"
    12  	"github.com/jingcheng-WU/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.EVCompute, a contains the
    23  // orthonormal eigenvectors of A on exit, otherwise jobz must be lapack.EVNone
    24  // and on exit the specified 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  	switch {
    32  	case jobz != lapack.EVNone && jobz != lapack.EVCompute:
    33  		panic(badEVJob)
    34  	case uplo != blas.Upper && uplo != blas.Lower:
    35  		panic(badUplo)
    36  	case n < 0:
    37  		panic(nLT0)
    38  	case lda < max(1, n):
    39  		panic(badLdA)
    40  	case lwork < max(1, 3*n-1) && lwork != -1:
    41  		panic(badLWork)
    42  	case len(work) < max(1, lwork):
    43  		panic(shortWork)
    44  	}
    45  
    46  	// Quick return if possible.
    47  	if n == 0 {
    48  		return true
    49  	}
    50  
    51  	var opts string
    52  	if uplo == blas.Upper {
    53  		opts = "U"
    54  	} else {
    55  		opts = "L"
    56  	}
    57  	nb := impl.Ilaenv(1, "DSYTRD", opts, n, -1, -1, -1)
    58  	lworkopt := max(1, (nb+2)*n)
    59  	if lwork == -1 {
    60  		work[0] = float64(lworkopt)
    61  		return
    62  	}
    63  
    64  	switch {
    65  	case len(a) < (n-1)*lda+n:
    66  		panic(shortA)
    67  	case len(w) < n:
    68  		panic(shortW)
    69  	}
    70  
    71  	if n == 1 {
    72  		w[0] = a[0]
    73  		work[0] = 2
    74  		if jobz == lapack.EVCompute {
    75  			a[0] = 1
    76  		}
    77  		return true
    78  	}
    79  
    80  	safmin := dlamchS
    81  	eps := dlamchP
    82  	smlnum := safmin / eps
    83  	bignum := 1 / smlnum
    84  	rmin := math.Sqrt(smlnum)
    85  	rmax := math.Sqrt(bignum)
    86  
    87  	// Scale matrix to allowable range, if necessary.
    88  	anrm := impl.Dlansy(lapack.MaxAbs, uplo, n, a, lda, work)
    89  	scaled := false
    90  	var sigma float64
    91  	if anrm > 0 && anrm < rmin {
    92  		scaled = true
    93  		sigma = rmin / anrm
    94  	} else if anrm > rmax {
    95  		scaled = true
    96  		sigma = rmax / anrm
    97  	}
    98  	if scaled {
    99  		kind := lapack.LowerTri
   100  		if uplo == blas.Upper {
   101  			kind = lapack.UpperTri
   102  		}
   103  		impl.Dlascl(kind, 0, 0, 1, sigma, n, n, a, lda)
   104  	}
   105  	var inde int
   106  	indtau := inde + n
   107  	indwork := indtau + n
   108  	llwork := lwork - indwork
   109  	impl.Dsytrd(uplo, n, a, lda, w, work[inde:], work[indtau:], work[indwork:], llwork)
   110  
   111  	// For eigenvalues only, call Dsterf. For eigenvectors, first call Dorgtr
   112  	// to generate the orthogonal matrix, then call Dsteqr.
   113  	if jobz == lapack.EVNone {
   114  		ok = impl.Dsterf(n, w, work[inde:])
   115  	} else {
   116  		impl.Dorgtr(uplo, n, a, lda, work[indtau:], work[indwork:], llwork)
   117  		ok = impl.Dsteqr(lapack.EVComp(jobz), n, w, work[inde:], a, lda, work[indtau:])
   118  	}
   119  	if !ok {
   120  		return false
   121  	}
   122  
   123  	// If the matrix was scaled, then rescale eigenvalues appropriately.
   124  	if scaled {
   125  		bi := blas64.Implementation()
   126  		bi.Dscal(n, 1/sigma, w, 1)
   127  	}
   128  	work[0] = float64(lworkopt)
   129  	return true
   130  }