diff --git a/lapack/gonum/dgeev.go b/lapack/gonum/dgeev.go index eb911d2a..9119a0d0 100644 --- a/lapack/gonum/dgeev.go +++ b/lapack/gonum/dgeev.go @@ -116,7 +116,7 @@ func (impl Implementation) Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob maxwrk := 2*n + n*impl.Ilaenv(1, "DGEHRD", " ", n, 1, n, 0) if wantvl || wantvr { maxwrk = max(maxwrk, 2*n+(n-1)*impl.Ilaenv(1, "DORGHR", " ", n, 1, n, -1)) - impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, 0, n-1, + impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, 0, n-1, nil, 1, nil, nil, nil, 1, work, -1) maxwrk = max(maxwrk, max(n+1, n+int(work[0]))) side := lapack.EVLeft @@ -176,7 +176,7 @@ func (impl Implementation) Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob impl.Dorghr(n, ilo, ihi, vl, ldvl, tau, work[iwrk:], lwork-iwrk) // Perform QR iteration, accumulating Schur vectors in VL. iwrk = n - first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi, + first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, ilo, ihi, a, lda, wr, wi, vl, ldvl, work[iwrk:], lwork-iwrk) if wantvr { // Want left and right eigenvectors. @@ -192,7 +192,7 @@ func (impl Implementation) Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob impl.Dorghr(n, ilo, ihi, vr, ldvr, tau, work[iwrk:], lwork-iwrk) // Perform QR iteration, accumulating Schur vectors in VR. iwrk = n - first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.OriginalEV, n, ilo, ihi, + first = impl.Dhseqr(lapack.EigenvaluesAndSchur, lapack.SchurOrig, n, ilo, ihi, a, lda, wr, wi, vr, ldvr, work[iwrk:], lwork-iwrk) } else { // Compute eigenvalues only. diff --git a/lapack/gonum/dhseqr.go b/lapack/gonum/dhseqr.go index 75c3b705..7e641b12 100644 --- a/lapack/gonum/dhseqr.go +++ b/lapack/gonum/dhseqr.go @@ -27,11 +27,11 @@ import ( // be computed. // For other values of job Dhseqr will panic. // -// If compz == lapack.None, no Schur vectors will be computed and Z will not be +// If compz == lapack.SchurNone, no Schur vectors will be computed and Z will not be // referenced. -// If compz == lapack.HessEV, on return Z will contain the matrix of Schur +// If compz == lapack.SchurHess, on return Z will contain the matrix of Schur // vectors of H. -// If compz == lapack.OriginalEV, on entry z is assumed to contain the orthogonal +// If compz == lapack.SchurOrig, on entry z is assumed to contain the orthogonal // matrix Q that is the identity except for the submatrix // Q[ilo:ihi+1,ilo:ihi+1]. On return z will be updated to the product Q*Z. // @@ -96,11 +96,11 @@ import ( // where U is an orthogonal matrix. The final H is upper Hessenberg and // H[unconverged:ihi+1,unconverged:ihi+1] is upper quasi-triangular. // -// If unconverged > 0 and compz == lapack.OriginalEV, then on return +// If unconverged > 0 and compz == lapack.SchurOrig, then on return // (final Z) = (initial Z) U, // where U is the orthogonal matrix in (*) regardless of the value of job. // -// If unconverged > 0 and compz == lapack.HessEV, then on return +// If unconverged > 0 and compz == lapack.SchurHess, then on return // (final Z) = U, // where U is the orthogonal matrix in (*) regardless of the value of job. // @@ -118,7 +118,7 @@ import ( // URL: http://dx.doi.org/10.1137/S0895479801384585 // // Dhseqr is an internal routine. It is exported for testing purposes. -func (impl Implementation) Dhseqr(job lapack.SchurJob, compz lapack.EVComp, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, z []float64, ldz int, work []float64, lwork int) (unconverged int) { +func (impl Implementation) Dhseqr(job lapack.SchurJob, compz lapack.SchurComp, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, z []float64, ldz int, work []float64, lwork int) (unconverged int) { var wantt bool switch job { default: @@ -130,9 +130,9 @@ func (impl Implementation) Dhseqr(job lapack.SchurJob, compz lapack.EVComp, n, i var wantz bool switch compz { default: - panic(badEVComp) - case lapack.None: - case lapack.HessEV, lapack.OriginalEV: + panic(badSchurComp) + case lapack.SchurNone: + case lapack.SchurHess, lapack.SchurOrig: wantz = true } switch { @@ -197,7 +197,7 @@ func (impl Implementation) Dhseqr(job lapack.SchurJob, compz lapack.EVComp, n, i } // Initialize Z to identity matrix if requested. - if compz == lapack.HessEV { + if compz == lapack.SchurHess { impl.Dlaset(blas.All, n, n, 0, 1, z, ldz) } diff --git a/lapack/gonum/dsteqr.go b/lapack/gonum/dsteqr.go index 0e1125e5..1da15a60 100644 --- a/lapack/gonum/dsteqr.go +++ b/lapack/gonum/dsteqr.go @@ -26,11 +26,11 @@ import ( // Dsteqr will panic otherwise. // // z, on entry, contains the n×n orthogonal matrix used in the reduction to -// tridiagonal form if compz == lapack.OriginalEV. On exit, if -// compz == lapack.OriginalEV, z contains the orthonormal eigenvectors of the -// original symmetric matrix, and if compz == lapack.TridiagEV, z contains the +// tridiagonal form if compz == lapack.EVOriginal. On exit, if +// compz == lapack.EVOriginal, z contains the orthonormal eigenvectors of the +// original symmetric matrix, and if compz == lapack.EVTridiag, z contains the // orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used -// if compz == lapack.None. +// if compz == lapack.EVCompNone. // // work must have length at least max(1, 2*n-2) if the eigenvectors are computed, // and Dsteqr will panic otherwise. @@ -46,10 +46,10 @@ func (impl Implementation) Dsteqr(compz lapack.EVComp, n int, d, e, z []float64, if len(e) < n-1 { panic(badE) } - if compz != lapack.None && compz != lapack.TridiagEV && compz != lapack.OriginalEV { + if compz != lapack.EVCompNone && compz != lapack.EVTridiag && compz != lapack.EVOrig { panic(badEVComp) } - if compz != lapack.None { + if compz != lapack.EVCompNone { if len(work) < max(1, 2*n-2) { panic(badWork) } @@ -57,9 +57,9 @@ func (impl Implementation) Dsteqr(compz lapack.EVComp, n int, d, e, z []float64, } var icompz int - if compz == lapack.OriginalEV { + if compz == lapack.EVOrig { icompz = 1 - } else if compz == lapack.TridiagEV { + } else if compz == lapack.EVTridiag { icompz = 2 } diff --git a/lapack/gonum/dtrexc.go b/lapack/gonum/dtrexc.go index 40e03785..1953fca9 100644 --- a/lapack/gonum/dtrexc.go +++ b/lapack/gonum/dtrexc.go @@ -18,8 +18,8 @@ import "gonum.org/v1/gonum/lapack" // as Z^T*T*Z, and will be again in Schur canonical form. // // If compq is lapack.UpdateSchur, on return the matrix Q of Schur vectors will be -// updated by postmultiplying it with Z. -// If compq is lapack.None, the matrix Q is not referenced and will not be +// updated by post-multiplying it with Z. +// If compq is lapack.UpdateSchurNone, the matrix Q is not referenced and will not be // updated. // For other values of compq Dtrexc will panic. // @@ -45,13 +45,13 @@ import "gonum.org/v1/gonum/lapack" // work must have length at least n, otherwise Dtrexc will panic. // // Dtrexc is an internal routine. It is exported for testing purposes. -func (impl Implementation) Dtrexc(compq lapack.EVComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool) { +func (impl Implementation) Dtrexc(compq lapack.UpdateSchurComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool) { checkMatrix(n, n, t, ldt) var wantq bool switch compq { default: panic("lapack: bad value of compq") - case lapack.None: + case lapack.UpdateSchurNone: // Nothing to do because wantq is already false. case lapack.UpdateSchur: wantq = true diff --git a/lapack/gonum/general.go b/lapack/gonum/general.go index 4615c52b..d8fb0c7d 100644 --- a/lapack/gonum/general.go +++ b/lapack/gonum/general.go @@ -45,6 +45,7 @@ const ( badNorm = "lapack: bad norm" badPivot = "lapack: bad pivot" badS = "lapack: s has insufficient length" + badSchurComp = "lapack: bad SchurComp" badSchurJob = "lapack: bad SchurJob" badShifts = "lapack: bad shifts" badSide = "lapack: bad side" diff --git a/lapack/lapack.go b/lapack/lapack.go index 45943bdf..cad3b9ca 100644 --- a/lapack/lapack.go +++ b/lapack/lapack.go @@ -127,23 +127,13 @@ const ( GSVDNone GSVDJob = 'N' // Do not compute orthogonal matrix ) -// EVComp specifies how eigenvectors are computed. +// EVComp specifies how eigenvectors are computed in Dsteqr. type EVComp byte const ( - // OriginalEV specifies to compute the eigenvectors of the original - // matrix. - OriginalEV EVComp = 'V' - // TridiagEV specifies to compute both the eigenvectors of the input - // tridiagonal matrix. - TridiagEV EVComp = 'I' - // HessEV specifies to compute both the eigenvectors of the input upper - // Hessenberg matrix. - HessEV EVComp = 'I' - - // UpdateSchur specifies that the matrix of Schur vectors will be - // updated by Dtrexc. - UpdateSchur EVComp = 'V' + EVOrig EVComp = 'V' // Compute eigenvectors of the original symmetric matrix. + EVTridiag EVComp = 'I' // Compute eigenvectors of the tridiagonal matrix. + EVCompNone EVComp = 'N' // Do not compute eigenvectors. ) // EVJob specifies whether eigenvectors are computed in Dsyev. @@ -188,6 +178,23 @@ const ( EigenvaluesAndSchur SchurJob = 'S' ) +// SchurComp specifies whether and how the Schur vectors are computed in Dhseqr. +type SchurComp byte + +const ( + SchurNone SchurComp = 'N' // Schur vectors are not computed. + SchurHess SchurComp = 'I' // Schur vectors of the upper Hessenberg marix are computed. + SchurOrig SchurComp = 'V' // Schur vectors of the original matrix are computed. +) + +// UpdateSchurComp specifies whether the matrix of Schur vectors is updated in Dtrexc. +type UpdateSchurComp byte + +const ( + UpdateSchur UpdateSchurComp = 'V' // The matrix of Schur vectors is updated. + UpdateSchurNone UpdateSchurComp = 'N' // The matrix of Schur vectors is not updated. +) + // EVSide specifies what eigenvectors are computed in Dtrevc3. type EVSide byte diff --git a/lapack/testlapack/dhseqr.go b/lapack/testlapack/dhseqr.go index b3866927..c7d0830e 100644 --- a/lapack/testlapack/dhseqr.go +++ b/lapack/testlapack/dhseqr.go @@ -16,7 +16,7 @@ import ( ) type Dhseqrer interface { - Dhseqr(job lapack.SchurJob, compz lapack.EVComp, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, + Dhseqr(job lapack.SchurJob, compz lapack.SchurComp, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, z []float64, ldz int, work []float64, lwork int) int } @@ -57,11 +57,11 @@ func testDhseqr(t *testing.T, impl Dhseqrer, i int, test dhseqrTest, job lapack. copyGeneral(h, blas64.General{Rows: n, Cols: n, Stride: max(1, n), Data: test.h}) hCopy := cloneGeneral(h) - var compz lapack.EVComp = lapack.None + compz := lapack.SchurNone z := blas64.General{Stride: max(1, n)} if wantz { // First, let Dhseqr initialize Z to the identity matrix. - compz = lapack.HessEV + compz = lapack.SchurHess z = nanGeneral(n, n, n+extra) } @@ -70,7 +70,7 @@ func testDhseqr(t *testing.T, impl Dhseqrer, i int, test dhseqrTest, job lapack. work := nanSlice(max(1, n)) if optwork { - impl.Dhseqr(job, lapack.HessEV, n, ilo, ihi, nil, h.Stride, nil, nil, nil, z.Stride, work, -1) + impl.Dhseqr(job, lapack.SchurHess, n, ilo, ihi, nil, h.Stride, nil, nil, nil, z.Stride, work, -1) work = nanSlice(int(work[0])) } @@ -196,7 +196,7 @@ func testDhseqr(t *testing.T, impl Dhseqrer, i int, test dhseqrTest, job lapack. copyGeneral(h, hCopy) // Call Dhseqr again with the identity matrix given explicitly in Q. q := eye(n, n+extra) - impl.Dhseqr(job, lapack.OriginalEV, n, ilo, ihi, h.Data, h.Stride, wr, wi, q.Data, q.Stride, work, len(work)) + impl.Dhseqr(job, lapack.SchurOrig, n, ilo, ihi, h.Data, h.Stride, wr, wi, q.Data, q.Stride, work, len(work)) if !equalApproxGeneral(z, q, 0) { t.Errorf("%v: Z and Q are not equal", prefix) } diff --git a/lapack/testlapack/dsteqr.go b/lapack/testlapack/dsteqr.go index 55978645..f95a7419 100644 --- a/lapack/testlapack/dsteqr.go +++ b/lapack/testlapack/dsteqr.go @@ -22,7 +22,7 @@ type Dsteqrer interface { func DsteqrTest(t *testing.T, impl Dsteqrer) { rnd := rand.New(rand.NewSource(1)) - for _, compz := range []lapack.EVComp{lapack.OriginalEV, lapack.TridiagEV} { + for _, compz := range []lapack.EVComp{lapack.EVOrig, lapack.EVTridiag} { for _, test := range []struct { n, lda int }{ @@ -59,7 +59,7 @@ func DsteqrTest(t *testing.T, impl Dsteqrer) { copy(eCopy, e) aCopy := make([]float64, len(a)) copy(aCopy, a) - if compz == lapack.OriginalEV { + if compz == lapack.EVOrig { uplo := blas.Upper tau := make([]float64, n) work := make([]float64, 1) @@ -92,7 +92,7 @@ func DsteqrTest(t *testing.T, impl Dsteqrer) { copy(dAns, d) var truth blas64.General - if compz == lapack.OriginalEV { + if compz == lapack.EVOrig { truth = blas64.General{ Rows: n, Cols: n, @@ -130,7 +130,7 @@ func DsteqrTest(t *testing.T, impl Dsteqrer) { } if !eigenDecompCorrect(d, truth, V) { t.Errorf("Eigen reconstruction mismatch. fromFull = %v, n = %v", - compz == lapack.OriginalEV, n) + compz == lapack.EVOrig, n) } // Compare eigenvalues when not computing eigenvectors. diff --git a/lapack/testlapack/dtrexc.go b/lapack/testlapack/dtrexc.go index dedc106c..f684f9ad 100644 --- a/lapack/testlapack/dtrexc.go +++ b/lapack/testlapack/dtrexc.go @@ -17,13 +17,13 @@ import ( ) type Dtrexcer interface { - Dtrexc(compq lapack.EVComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool) + Dtrexc(compq lapack.UpdateSchurComp, n int, t []float64, ldt int, q []float64, ldq int, ifst, ilst int, work []float64) (ifstOut, ilstOut int, ok bool) } func DtrexcTest(t *testing.T, impl Dtrexcer) { rnd := rand.New(rand.NewSource(1)) - for _, compq := range []lapack.EVComp{lapack.None, lapack.UpdateSchur} { + for _, compq := range []lapack.UpdateSchurComp{lapack.UpdateSchurNone, lapack.UpdateSchur} { for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 53} { for _, extra := range []int{0, 1, 11} { for cas := 0; cas < 100; cas++ { @@ -36,7 +36,7 @@ func DtrexcTest(t *testing.T, impl Dtrexcer) { } } - for _, compq := range []lapack.EVComp{lapack.None, lapack.UpdateSchur} { + for _, compq := range []lapack.UpdateSchurComp{lapack.UpdateSchurNone, lapack.UpdateSchur} { for _, extra := range []int{0, 1, 11} { tmat := randomSchurCanonical(0, extra, rnd) testDtrexc(t, impl, compq, tmat, 0, 0, extra, rnd) @@ -44,7 +44,7 @@ func DtrexcTest(t *testing.T, impl Dtrexcer) { } } -func testDtrexc(t *testing.T, impl Dtrexcer, compq lapack.EVComp, tmat blas64.General, ifst, ilst, extra int, rnd *rand.Rand) { +func testDtrexc(t *testing.T, impl Dtrexcer, compq lapack.UpdateSchurComp, tmat blas64.General, ifst, ilst, extra int, rnd *rand.Rand) { const tol = 1e-13 n := tmat.Rows