github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/gonum/blas"
    11  	"github.com/gonum/blas/blas64"
    12  	"github.com/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  }