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
5package gonum
6
7import (
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// Dgeev computes the eigenvalues and, optionally, the left and/or right
16// eigenvectors for an n×n real nonsymmetric matrix A.
17//
18// The right eigenvector v_j of A corresponding to an eigenvalue λ_j
19// is defined by
20//  A v_j = λ_j v_j,
21// and the left eigenvector u_j corresponding to an eigenvalue λ_j is defined by
22//  u_jᴴ A = λ_j u_jᴴ,
23// where u_jᴴ is the conjugate transpose of u_j.
24//
25// On return, A will be overwritten and the left and right eigenvectors will be
26// stored, respectively, in the columns of the n×n matrices VL and VR in the
27// same order as their eigenvalues. If the j-th eigenvalue is real, then
28//  u_j = VL[:,j],
29//  v_j = VR[:,j],
30// and if it is not real, then j and j+1 form a complex conjugate pair and the
31// eigenvectors can be recovered as
32//  u_j     = VL[:,j] + i*VL[:,j+1],
33//  u_{j+1} = VL[:,j] - i*VL[:,j+1],
34//  v_j     = VR[:,j] + i*VR[:,j+1],
35//  v_{j+1} = VR[:,j] - i*VR[:,j+1],
36// where i is the imaginary unit. The computed eigenvectors are normalized to
37// have Euclidean norm equal to 1 and largest component real.
38//
39// Left eigenvectors will be computed only if jobvl == lapack.LeftEVCompute,
40// otherwise jobvl must be lapack.LeftEVNone.
41// Right eigenvectors will be computed only if jobvr == lapack.RightEVCompute,
42// otherwise jobvr must be lapack.RightEVNone.
43// For other values of jobvl and jobvr Dgeev will panic.
44//
45// wr and wi contain the real and imaginary parts, respectively, of the computed
46// eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with
47// the eigenvalue having the positive imaginary part first.
48// wr and wi must have length n, and Dgeev will panic otherwise.
49//
50// work must have length at least lwork and lwork must be at least max(1,4*n) if
51// the left or right eigenvectors are computed, and at least max(1,3*n) if no
52// eigenvectors are computed. For good performance, lwork must generally be
53// larger.  On return, optimal value of lwork will be stored in work[0].
54//
55// If lwork == -1, instead of performing Dgeev, the function only calculates the
56// optimal vaule of lwork and stores it into work[0].
57//
58// On return, first is the index of the first valid eigenvalue. If first == 0,
59// all eigenvalues and eigenvectors have been computed. If first is positive,
60// Dgeev failed to compute all the eigenvalues, no eigenvectors have been
61// computed and wr[first:] and wi[first:] contain those eigenvalues which have
62// converged.
63func (impl Implementation) Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, n int, a []float64, lda int, wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) (first int) {
64	wantvl := jobvl == lapack.LeftEVCompute
65	wantvr := jobvr == lapack.RightEVCompute
66	var minwrk int
67	if wantvl || wantvr {
68		minwrk = max(1, 4*n)
69	} else {
70		minwrk = max(1, 3*n)
71	}
72	switch {
73	case jobvl != lapack.LeftEVCompute && jobvl != lapack.LeftEVNone:
74		panic(badLeftEVJob)
75	case jobvr != lapack.RightEVCompute && jobvr != lapack.RightEVNone:
76		panic(badRightEVJob)
77	case n < 0:
78		panic(nLT0)
79	case lda < max(1, n):
80		panic(badLdA)
81	case ldvl < 1 || (ldvl < n && wantvl):
82		panic(badLdVL)
83	case ldvr < 1 || (ldvr < n && wantvr):
84		panic(badLdVR)
85	case lwork < minwrk && lwork != -1:
86		panic(badLWork)
87	case len(work) < lwork:
88		panic(shortWork)
89	}
90
91	// Quick return if possible.
92	if n == 0 {
93		work[0] = 1
94		return 0
95	}
96
97	maxwrk := 2*n + n*impl.Ilaenv(1, "DGEHRD", " ", n, 1, n, 0)
98	if wantvl || wantvr {
99		maxwrk = max(maxwrk, 2*n+(n-1)*impl.Ilaenv(1, "DORGHR", " ", n, 1, n, -1))
100		impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, 0, n-1,
101			a, lda, wr, wi, nil, n, work, -1)
102		maxwrk = max(maxwrk, max(n+1, n+int(work[0])))
103		side := lapack.EVLeft
104		if wantvr {
105			side = lapack.EVRight
106		}
107		impl.Dtrevc3(side, lapack.EVAllMulQ, nil, n, a, lda, vl, ldvl, vr, ldvr,
108			n, work, -1)
109		maxwrk = max(maxwrk, n+int(work[0]))
110		maxwrk = max(maxwrk, 4*n)
111	} else {
112		impl.Dhseqr(lapack.EigenvaluesOnly, lapack.SchurNone, n, 0, n-1,
113			a, lda, wr, wi, vr, ldvr, work, -1)
114		maxwrk = max(maxwrk, max(n+1, n+int(work[0])))
115	}
116	maxwrk = max(maxwrk, minwrk)
117
118	if lwork == -1 {
119		work[0] = float64(maxwrk)
120		return 0
121	}
122
123	switch {
124	case len(a) < (n-1)*lda+n:
125		panic(shortA)
126	case len(wr) != n:
127		panic(badLenWr)
128	case len(wi) != n:
129		panic(badLenWi)
130	case len(vl) < (n-1)*ldvl+n && wantvl:
131		panic(shortVL)
132	case len(vr) < (n-1)*ldvr+n && wantvr:
133		panic(shortVR)
134	}
135
136	// Get machine constants.
137	smlnum := math.Sqrt(dlamchS) / dlamchP
138	bignum := 1 / smlnum
139
140	// Scale A if max element outside range [smlnum,bignum].
141	anrm := impl.Dlange(lapack.MaxAbs, n, n, a, lda, nil)
142	var scalea bool
143	var cscale float64
144	if 0 < anrm && anrm < smlnum {
145		scalea = true
146		cscale = smlnum
147	} else if anrm > bignum {
148		scalea = true
149		cscale = bignum
150	}
151	if scalea {
152		impl.Dlascl(lapack.General, 0, 0, anrm, cscale, n, n, a, lda)
153	}
154
155	// Balance the matrix.
156	workbal := work[:n]
157	ilo, ihi := impl.Dgebal(lapack.PermuteScale, n, a, lda, workbal)
158
159	// Reduce to upper Hessenberg form.
160	iwrk := 2 * n
161	tau := work[n : iwrk-1]
162	impl.Dgehrd(n, ilo, ihi, a, lda, tau, work[iwrk:], lwork-iwrk)
163
164	var side lapack.EVSide
165	if wantvl {
166		side = lapack.EVLeft
167		// Copy Householder vectors to VL.
168		impl.Dlacpy(blas.Lower, n, n, a, lda, vl, ldvl)
169		// Generate orthogonal matrix in VL.
170		impl.Dorghr(n, ilo, ihi, vl, ldvl, tau, work[iwrk:], lwork-iwrk)
171		// Perform QR iteration, accumulating Schur vectors in VL.
172		iwrk = n
173		first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, ilo, ihi,
174			a, lda, wr, wi, vl, ldvl, work[iwrk:], lwork-iwrk)
175		if wantvr {
176			// Want left and right eigenvectors.
177			// Copy Schur vectors to VR.
178			side = lapack.EVBoth
179			impl.Dlacpy(blas.All, n, n, vl, ldvl, vr, ldvr)
180		}
181	} else if wantvr {
182		side = lapack.EVRight
183		// Copy Householder vectors to VR.
184		impl.Dlacpy(blas.Lower, n, n, a, lda, vr, ldvr)
185		// Generate orthogonal matrix in VR.
186		impl.Dorghr(n, ilo, ihi, vr, ldvr, tau, work[iwrk:], lwork-iwrk)
187		// Perform QR iteration, accumulating Schur vectors in VR.
188		iwrk = n
189		first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, ilo, ihi,
190			a, lda, wr, wi, vr, ldvr, work[iwrk:], lwork-iwrk)
191	} else {
192		// Compute eigenvalues only.
193		iwrk = n
194		first = impl.Dhseqr(lapack.EigenvaluesOnly, lapack.SchurNone, n, ilo, ihi,
195			a, lda, wr, wi, nil, 1, work[iwrk:], lwork-iwrk)
196	}
197
198	if first > 0 {
199		if scalea {
200			// Undo scaling.
201			impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1)
202			impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1)
203			impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wr, 1)
204			impl.Dlascl(lapack.General, 0, 0, cscale, anrm, ilo, 1, wi, 1)
205		}
206		work[0] = float64(maxwrk)
207		return first
208	}
209
210	if wantvl || wantvr {
211		// Compute left and/or right eigenvectors.
212		impl.Dtrevc3(side, lapack.EVAllMulQ, nil, n,
213			a, lda, vl, ldvl, vr, ldvr, n, work[iwrk:], lwork-iwrk)
214	}
215	bi := blas64.Implementation()
216	if wantvl {
217		// Undo balancing of left eigenvectors.
218		impl.Dgebak(lapack.PermuteScale, lapack.EVLeft, n, ilo, ihi, workbal, n, vl, ldvl)
219		// Normalize left eigenvectors and make largest component real.
220		for i, wii := range wi {
221			if wii < 0 {
222				continue
223			}
224			if wii == 0 {
225				scl := 1 / bi.Dnrm2(n, vl[i:], ldvl)
226				bi.Dscal(n, scl, vl[i:], ldvl)
227				continue
228			}
229			scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vl[i:], ldvl), bi.Dnrm2(n, vl[i+1:], ldvl))
230			bi.Dscal(n, scl, vl[i:], ldvl)
231			bi.Dscal(n, scl, vl[i+1:], ldvl)
232			for k := 0; k < n; k++ {
233				vi := vl[k*ldvl+i]
234				vi1 := vl[k*ldvl+i+1]
235				work[iwrk+k] = vi*vi + vi1*vi1
236			}
237			k := bi.Idamax(n, work[iwrk:iwrk+n], 1)
238			cs, sn, _ := impl.Dlartg(vl[k*ldvl+i], vl[k*ldvl+i+1])
239			bi.Drot(n, vl[i:], ldvl, vl[i+1:], ldvl, cs, sn)
240			vl[k*ldvl+i+1] = 0
241		}
242	}
243	if wantvr {
244		// Undo balancing of right eigenvectors.
245		impl.Dgebak(lapack.PermuteScale, lapack.EVRight, n, ilo, ihi, workbal, n, vr, ldvr)
246		// Normalize right eigenvectors and make largest component real.
247		for i, wii := range wi {
248			if wii < 0 {
249				continue
250			}
251			if wii == 0 {
252				scl := 1 / bi.Dnrm2(n, vr[i:], ldvr)
253				bi.Dscal(n, scl, vr[i:], ldvr)
254				continue
255			}
256			scl := 1 / impl.Dlapy2(bi.Dnrm2(n, vr[i:], ldvr), bi.Dnrm2(n, vr[i+1:], ldvr))
257			bi.Dscal(n, scl, vr[i:], ldvr)
258			bi.Dscal(n, scl, vr[i+1:], ldvr)
259			for k := 0; k < n; k++ {
260				vi := vr[k*ldvr+i]
261				vi1 := vr[k*ldvr+i+1]
262				work[iwrk+k] = vi*vi + vi1*vi1
263			}
264			k := bi.Idamax(n, work[iwrk:iwrk+n], 1)
265			cs, sn, _ := impl.Dlartg(vr[k*ldvr+i], vr[k*ldvr+i+1])
266			bi.Drot(n, vr[i:], ldvr, vr[i+1:], ldvr, cs, sn)
267			vr[k*ldvr+i+1] = 0
268		}
269	}
270
271	if scalea {
272		// Undo scaling.
273		impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wr[first:], 1)
274		impl.Dlascl(lapack.General, 0, 0, cscale, anrm, n-first, 1, wi[first:], 1)
275	}
276
277	work[0] = float64(maxwrk)
278	return first
279}
280