Add wrapperfor ?pteqr

[mdhaber]

See NumFOCUS Small Development Grant Proposal 2019 Round 3 proposal for background.

Add wrapper for ?pteqr.

Suggested signature:

! d, e, z, work, info  = pteqr(d, e, z, compz="N", overwrite_d=0, overwrite_e=0, overwrite_z)

Test idea (compz="N"):

• Generate random tridiagonal SPD tridiagonal A. The fact that "a Hermitian strictly diagonally dominant matrix with real positive diagonal entries is positive definite" can be used to generate one. For instance, in the real case, if all the diagonal elements are >2 and the off-diagonals are <1, the matrix will be PD.
• Decompose A with ?pteqr
• Compare eigenvalues with those computed by scipy.linal.eigh

Test idea (compz="I"):

• Generate random tridiagonal SPD tridiagonal A.
• Decompose A with ?pteqr
• Compare eigenvalues with those computed by scipy.linal.eigh
• Verify that Z is orthogonal (Z@Z.T is identity matrix)
• Verify that Z @ diag(d) @ Z.T is the original matrix A.

I'll come back to tests for compz="V" later. Update: see #11224 (comment).

Also test for singular matrix, non-spd matrix, and incorrect/incompatible array sizes.

Also: implement all examples from the NAG manual as additional tests.

[mdhaber]

@swallan

[swallan]

Here is the branch that I started for the ?pteqr routine.

Something that I encountered when looking at the NAG examples that I questioned was that the NAG page said both d and e are Real (Kind=nag_wp) arrays, but it looks like the complex example here https://www.nag.com/numeric/fl/nagdoc_latest/html/f08/f08juf.html has some complex entries on e. The lapack page just says double precision array.

So, should I keep the randomly generated e real or allow it to be complex?

[mdhaber]

Allow it to be complex if it works.

[ilayn]

Try to follow Netlib documentation. NAG has no obligation to be LAPACK compatible and sometimes they have different flavors. That's why they have their own naming convention. Also NAG examples are low precision so if possible make random examples and use the NAG ones for validation.

[mdhaber]

@swallan I think you are you referring to the complex values in the example matrix A. e is always real.

Before calling zpteqr, the NAG example uses zhetrd to reduce the matrix to tridiagonal form, then zungtr forms an orthogonal matrix from factors returned by zhetrd. I think it then passes a real tridiagonal matrix and complex orthogonal matrix into zpteqr. So please pass only real d and e into the routines, even for complex data types.

I wrote above:

I'll come back to tests for compz="V" later.

That's where the complex data types would actually be used. Based on the NAG examples and documentation, see if you can figure out how tests for compz="V" would work. In short, you will start with a full (not tridiagonal) symmetric/Hessian matrix A, convert it to tridiagonal form (with ?sytrd and ?orgtr or their complex equivalents), and then use ?pteqr. You'd then compare the eigenvalues of the original matrix A produced by ?pteqr and scipy.linalg.eigh. If you have trouble, we can work on it in office hours.

[swallan]

@mdhaber Yes, I was referring the complex A– I will take a crack at it and see where I can get.

[Kai-Striega]

@swallan I've pushed the wrapper to my branch lapack_support_pteqr. I know that I have been a bit behind over the last week(s), but I'm going to be able to get these out fairly quickly from now.

@mdhaber said that you have been working on these on your own branch, do you have any which are ready/almost ready to go? I can prioritise those, so we can start getting these merged.

[swallan]

I'm not too picky–they're all mostly written.

[swallan]

@Kai-Striega
I'm taking a look at running the tests for this one, and I think that in the wrapper z cannot depend on n or should accept None since z is not used in all of the compz values -

If COMPZ = 'N', then Z is not referenced.

And my reading of the documentation makes me think that z is not an input when compz=I either

[Kai-Striega]

I think it's still included, it's just not referenced. For the case where compz=I, z is definitely included as, on return, it contains the orthonormal eigenvectors of the tridiagonal matrix. Although I agree that not having to explicitly pass z would be nice.

…

On Mon., 24 Feb. 2020, 13:34 Sam Wallan, ***@***.***> wrote: @Kai-Striega <https://github.com/Kai-Striega> I'm taking a look at running the tests for this one, and I think that in the wrapper z cannot depend on n or should accept None since z is not used in all of the compz values - If COMPZ = 'N', then Z is not referenced. And my reading of the documentation makes me think that z is not an input when compz=I either — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11224?email_source=notifications&email_token=AHFOJNLMXCG7JPREHG7WTXDRENL7NA5CNFSM4J3CLJM2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMWVJ5Q#issuecomment-590173430>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHFOJNIXTUNSWAWCTI6SJQDRENL7NANCNFSM4J3CLJMQ> .

[swallan]

Opened a pull request to add the tests for this wrapper.

[Kai-Striega]

@ilayn I'd like to check something with you. The definition of ldz is:

LDZ is INTEGER
The leading dimension of the array Z. LDZ >= 1, and if COMPZ = 'V' or 'I', LDZ >= max(1,N).

This means that ldz can be either 1 or n, depending on compz. One way to set this up is:

integer intent(hide), depend(z, n, compz), check(ldz>=((*compz=='N')?1:n)) :: ldz = max(1, shape(z, 0))

OR, moving the check into the definition:

integer intent(hide), depend(z, n, compz) :: ldz = max((*compz=='N')?1:n, shape(z, 0))

Which works... but looks quite convoluted and error prone. Is there a better way to achieve the same result?

[ilayn]

Ideally if the char is actually acting like a switch we convert it into integers 0 1 2 to make the switch easier and comparisons to be over integers instead of chars. You can find many examples of compute_v, compute_z etc.

Here the switch is a boolean like. If compute_z is 0 z is an empty array hence optional. But if not zero it needs a check of being n x n. Same for ldz if 0 ldz=1 or it is n.

[Kai-Striega]

Thanks I'll take a look at compute_z etc. tomorrow morning!