ENH: Wrap LAPACK's dsytrd #7780

nschloe · 2017-08-22T08:27:38Z

This PR adds the wrapper for LAPACK's dsytrd, a routine that takes a real symmetric matrix A and transforms it to a tridiagonal matrix by an orthogonal similarity tranformation

Q^T * A * Q = T

With test.

I would have wrapped its complex variant *HETRD as well, but I didn't know how to map the complex type <ctype2c> to its corresponding real type in the wrapper code. (HETRD transforms a complex Hermitian matrix to a real symmetric tridiagonal matrix.)

Fixes #7775.

ilayn

I think once the <prefix> part is fixed this can build and the tests can be evaluated properly. One thing to look out for is not to cause segfaults even when used inappropriately. Recently we had a case with supplying np.empty((0,0)) was causing segfaults so better add it and possibly others if any, to the tests.

ilayn · 2017-08-22T09:56:33Z

scipy/linalg/flapack.pyf.src

+     integer intent(in) :: n
+     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+
+     integer intent(hide),depend(n) :: lda=n


Haven't checked if lwork=-1 causes a shortcut in the Fortran source, but still better to keep lda = MAX(n,1) to be on the safe side as it causes trouble on the Fortran side with zero-sized arrays. n=0 OK for LAPACK for which they have often a safeguard as a convention but lda=0 is not.

ilayn · 2017-08-22T09:57:31Z

scipy/linalg/flapack.pyf.src

+     callstatement (*f2py_func)((lower?"L":"U"),&n,a,&lda,d,e,tau,work,&lwork,&info);
+     callprotoargument char*,int*,<ctype2>*,int*,<ctype2>*,<ctype2>*,<ctype2>*,<ctype2>*,int*,int*
+
+     integer intent(hide),depend(a) :: lda=shape(a,0)


Use lda = MAX(shape(a,0),1) to avoid lda=0

ilayn · 2017-08-22T10:03:10Z

scipy/linalg/tests/test_lapack.py

+        A = np.random.normal(size=(n, n))
+        A += A.T
+
+        sytrd, sytrd_lwork = get_lapack_funcs(('sytrd', 'sytrd_lwork'))


You would add the array you want to get the LAPACK routine for as get_lapack_funcs(('sytrd', 'sytrd_lwork'), (A,)) such that d is selected. Same should also test for s flavor with a float32 array.

ilayn · 2017-08-22T10:08:18Z

scipy/linalg/flapack.pyf.src

+     <ftype2> intent(out) :: work
+     integer intent(hide) :: lwork = -1
+     integer intent(out) :: info
+   end subroutine <prefix>sytrd_lwork


This needs to be <prefix2> instead.

ilayn · 2017-08-22T10:13:51Z

scipy/linalg/tests/test_lapack.py

+
+        # query lwork
+        lwork, info = sytrd_lwork(n)
+        assert info == 0


Naked asserts are not allowed. You can use assert_equal from numpy.testing module

ilayn · 2017-08-22T10:14:07Z

scipy/linalg/tests/test_lapack.py

+        assert info == 0
+
+        data, d, e, tau, info = sytrd(A, lwork)
+        assert info == 0


Same comment about assert

pv · 2017-08-22T12:34:30Z

scipy/linalg/tests/test_lapack.py

+        np.random.seed(42)
+        n = 3
+        A = np.random.normal(size=(n, n))
+        A += A.T


Inplace addition of overlapping arrays is undefined in Numpy < 1.13

nschloe · 2017-08-22T13:46:35Z

@ilayn When supplying a np.empty((0, 0)) matrix, the method fails with

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,)

that sounds pretty reasonable, right?

nschloe · 2017-08-23T08:57:51Z

I'm going to need some help here. The error message is

E           AssertionError: 
E           Not equal to tolerance rtol=0.000119209, atol=0
E           
E           (mismatch 22.22222222222223%)
E            x: array([[  2.941178e-02,   3.823532e-01,  -2.384186e-07],
E                  [  3.823531e-01,   4.970589e+00,  -5.830952e+00],
E                  [ -2.384186e-07,  -5.830952e+00,   6.000000e+00]], dtype=float32)
E            y: array([[ 0.029412,  0.382353,  0.      ],
E                  [ 0.382353,  4.970587, -5.830952],
E                  [ 0.      , -5.830952,  6.      ]], dtype=float32)

and really x and y look equal. Why is there a mismatch?

ev-br · 2017-08-23T09:14:08Z

What you can do is to drop into a pdb session (insert import pdb; pdb.set_trace() just before the assertion and print the difference

nschloe · 2017-08-23T09:22:37Z

@ev-br How to I execute just the LAPACK tests locally? (Otherwise it always takes half an hour.)

ev-br · 2017-08-23T09:36:28Z

$ python runtests.py -s linalg or $ python runtests.py -t /path/to/test/file

nschloe · 2017-08-24T14:43:19Z

@everyone Tests are passing now. Anything else for me to do?

ilayn

I've added some minor comments. I think when they are addressed, this is good to go.

ilayn · 2017-08-25T10:18:52Z

scipy/linalg/flapack.pyf.src

+     integer intent(hide),depend(a) :: n=shape(a,1)
+     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+
+     <ftype2> dimension(lda,n),intent(in,out,copy,out=c) :: a


This needs a squareness check before LAPACK chews it so it is best if we intervene

<ftype2> dimension(lda,n),check(shape(a,0)==shape(a,1)),intent(in,out,copy,out=c) :: a

ilayn · 2017-08-25T10:22:37Z

scipy/linalg/flapack.pyf.src

+     <ftype2> dimension(n),intent(out),depend(n) :: d
+     <ftype2> dimension(n-1),intent(out),depend(n) :: e
+     <ftype2> dimension(n-1),intent(out),depend(n) :: tau
+     integer intent(in),optional,depend(n),check(lwork>=MAX(n,1)) :: lwork = MAX(n,1)


Lwork is not needed to be n or higher. It just needs to be greater than or equal to 1. For this you can keep the default but the check needs to be as

integer intent(in),optional,depend(n),check(lwork>=1||lwork==-1) :: lwork = MAX(n,1)

Here we already define sytrd_lwork for -1 query, but some people do it manually so better to keep it general.

ilayn · 2017-08-25T10:41:12Z

scipy/linalg/tests/test_lapack.py

+    def test_sytrd(self):
+        for dtype in REAL_DTYPES:
+            # "random" symmetric matrix
+            A = np.array([


If you create A via

n = 3 A = np.zeros((n, n), dtype=dtype) A[np.triu_indices_from(A)] = np.arange(1, 2*n+1, dtype=dtype)

(or typing out manually if you wish), you can quickly test for lower=1 behavior for n of your choice since it will give the matrix and it's diagonal back together with zero arrays for a diagonal matrix. Then it can continue with lower=0 case for the precise calculation.

The comment can be worded as # some upper triangular array to avoid confusion about the randomness too.

ilayn · 2017-08-25T10:43:34Z

scipy/linalg/tests/test_lapack.py

+            # disable rtol here since some values in QTAQ and T are very close
+            # to 0.
+            assert_allclose(QTAQ, T, atol=5*np.finfo(dtype).eps, rtol=1.0)
+        return


I think we don't need this return here, no?

ilayn · 2017-08-25T17:23:55Z

I have just seen your final remark in the original PR comment. Let's not waste this opportunity to wrap the complex versions. This should do the (C/Z)HETRD versions

   subroutine <prefix2c>hetrd(lower,n,a,lda,d,e,tau,work,lwork,info)
     ! Reduce a complex hermitian matrix A to real symmetric
     ! tridiagonal form T by an orthogonal similarity transformation:
     ! Q**H * A * Q = T.

     callstatement (*f2py_func)((lower?"L":"U"),&n,a,&lda,d,e,tau,work,&lwork,&info);
     callprotoargument char*,int*,<ctype2c>*,int*,<ctype2>*,<ctype2>*,<ctype2c>*,<ctype2c>*,int*,int*

     integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
     integer intent(hide),depend(a) :: n=shape(a,1)
     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0

     <ftype2c> dimension(lda,n),check(shape(a,0)==shape(a,1)),intent(in,out,copy,out=c) :: a
     <ftype2> dimension(n),intent(out),depend(n) :: d
     <ftype2> dimension(n-1),intent(out),depend(n) :: e
     <ftype2c> dimension(n-1),intent(out),depend(n) :: tau
     integer intent(in),optional,depend(n),check(lwork>=1||lwork==-1) :: lwork = MAX(n,1)
     <ftype2c> dimension(lwork),intent(cache,hide) :: work
     integer intent(out) :: info
   end subroutine <prefix2c>hetrd

   subroutine <prefix2c>hetrd_lwork(lower,n,a,lda,d,e,tau,work,lwork,info)
     ! lwork computation for hetrd
     fortranname <prefix2c>hetrd
     callstatement (*f2py_func)((lower?"L":"U"),&n,&a,&lda,&d,&e,&tau,&work,&lwork,&info);
     callprotoargument char*,int*,<ctype2c>*,int*,<ctype2>*,<ctype2>*,<ctype2c>*,<ctype2c>*,int*,int*

     integer intent(in) :: n
     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0

     integer intent(hide),depend(n) :: lda=MAX(n,1)
     <ftype2c> intent(hide) :: a
     <ftype2> intent(hide) :: d
     <ftype2> intent(hide) :: e
     <ftype2c> intent(hide) :: tau
     <ftype2c> intent(out) :: work
     integer intent(hide) :: lwork = -1
     integer intent(out) :: info
   end subroutine <prefix2c>hetrd_lwork

together with their definitions in the lapack.py file

chetrd
zhetrd

chetrd_lwork
zhetrd_lwork

It would be great if you can extend the test to cover these too while we are at it.

eric-wieser · 2017-08-27T12:46:41Z

What does this do on (0,0) arrays? Does it still error?

Having fixed this type of bug in numpy, the error is typically that it should be lda = max(n, 1) as it says in the docs, not lda = n.

ilayn · 2017-08-27T12:58:33Z

@eric-wieser Yes it produces an error (also for 1x1 array). Since f2py needs to allocate arrays at the outset, the subdiagonals needs to be provided and in these cases subdiagonal sizes are -1 and 0 respectively which is problematic in both cases so it trips at the f2py check() condition with

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,)

In the other routines, we have also fixed a similar bug recently with the same conclusion (#7633)

eric-wieser · 2017-08-27T13:15:25Z

Doesn't seem like (1,1) should be a problem to me - isn't (0,) a valid size? I see your point about (0,0) though.

Either way, it might be nice to add tests for these cases.

nschloe · 2017-08-27T18:22:18Z

Doesn't seem like (1,1) should be a problem to me

I agree. I'll check it out tomorrow.

nschloe · 2017-08-28T07:34:39Z

The error that is produced when dealing with a 1x1 array is

>           data, d, e, tau, info = sytrd(A, lower=1)
E           ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)

I don't understand what variable gets (0,) – everything is protected with MAX(n, 1). Any ideas?

nschloe · 2017-08-28T07:59:31Z

Okay, so it's apparently about the arrays

     <ftype2> dimension(n-1),intent(out),depend(n) :: e
     <ftype2> dimension(n-1),intent(out),depend(n) :: tau

(representing the super/subdiagonal, e, and a list of Householder coefficients, tau). If n==1, it tries to allocate arrays of dimension 0, which as far as I know is alright with Fortran. SciPy's preprocessor seems to disagree though.

nschloe · 2017-08-28T08:31:07Z

Alright, so this appears to be a NumPy bug (numpy/numpy#9617). It wouldn't make much sense to add a test for n==1 at this point since the behavior will soon be corrected (I'm working on a NumPy PR right now). I've added a test case for n==0 to make sure that a ValueError is raised.

eric-wieser · 2017-08-28T11:10:38Z

I'd maybe argue for a dec.knownFailure test, or at least a commented out assertion with a TODO (referencing the numpy issue), so that testing it isn't completely forgotten.

nschloe · 2017-08-31T08:03:23Z

I've added a TODO note (although I'd consider this an upstream issue).

ilayn · 2017-08-31T08:30:58Z

Alright, so this appears to be a NumPy bug (numpy/numpy#9617). It wouldn't make much sense to add a test for n==1 at this point since the behavior will soon be corrected

Is it going to be backported?

nschloe · 2017-08-31T08:39:45Z

Is it going to be backported?

Probably a question to ask on the numpy side. Either way, no change is required on scipy for this to work properly.

pv · 2017-09-03T06:28:49Z

thanks, lgtm. I don't see a nice way to support the 1x1 case before f2py supports size-0 arrays.

nschloe added 5 commits August 21, 2017 15:58

add interface for LAPACK's sytrd

7cfe12a

fix sytrd interface, add test

153f313

add lwork computation for sytrd

1344428

simplify sytrd a bit

4fc1c07

better sytrd_lwork

c630268

nschloe mentioned this pull request Aug 22, 2017

wrap DSYTRD #7775

Closed

ilayn reviewed Aug 22, 2017

View reviewed changes

pv reviewed Aug 22, 2017

View reviewed changes

sytrd: address comments from code review

27e8a8c

nschloe added 4 commits August 22, 2017 15:46

sytrd: fix lwork specification

d7485a8

fix sytrd test

3db5b99

relax sytrd test tolerance

9ec8707

relax sytrd test tolerance again

ff4a8d0

sytrd: fix test tolerance

703033d

rgommers added enhancement A new feature or improvement scipy.linalg labels Aug 24, 2017

ilayn requested changes Aug 25, 2017

View reviewed changes

nschloe added 4 commits August 26, 2017 23:53

some small corrections for sytrd

4dbf36e

better sytrd tests

9bfef3a

add LAPACK wrappers for {c,z}hetrd

15aac40

add test for {c,z}hetrd

189104d

ilayn approved these changes Aug 27, 2017

View reviewed changes

{sy,he}trd tests: better array initialization

dae74bb

{sy,he}trd: missing depend(lwork)

993b45c

This was referenced Aug 28, 2017

Fortran preprocessor falsely disallows arrays of dimension 0 #7801

Closed

f2py: Fortran preprocessor falsely disallows arrays of dimension 0 numpy/numpy#9617

Closed

{sy,he}trd: assert_raises for 0x0 arrays

600c6ed

{sy,he}trd: add TODO note on tests with n=1

dc25c90

nschloe changed the title ~~Wrap LAPACK's dsytrd~~ ENH: Wrap LAPACK's dsytrd Aug 31, 2017

pv merged commit 892a9c4 into scipy:master Sep 3, 2017

pv added this to the 1.0.0 milestone Sep 3, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: Wrap LAPACK's dsytrd #7780

ENH: Wrap LAPACK's dsytrd #7780

nschloe commented Aug 22, 2017

ilayn left a comment

ilayn Aug 22, 2017

ilayn Aug 22, 2017

ilayn Aug 22, 2017

ilayn Aug 22, 2017

ilayn Aug 22, 2017

ilayn Aug 22, 2017

pv Aug 22, 2017

nschloe commented Aug 22, 2017

nschloe commented Aug 23, 2017

ev-br commented Aug 23, 2017

nschloe commented Aug 23, 2017

ev-br commented Aug 23, 2017

nschloe commented Aug 24, 2017

ilayn left a comment

ilayn Aug 25, 2017

ilayn Aug 25, 2017

ilayn Aug 25, 2017

ilayn Aug 25, 2017

ilayn Aug 25, 2017

ilayn commented Aug 25, 2017

eric-wieser commented Aug 27, 2017

ilayn commented Aug 27, 2017

eric-wieser commented Aug 27, 2017

nschloe commented Aug 27, 2017

nschloe commented Aug 28, 2017

nschloe commented Aug 28, 2017 •

edited

Loading

nschloe commented Aug 28, 2017

eric-wieser commented Aug 28, 2017

nschloe commented Aug 31, 2017

ilayn commented Aug 31, 2017

nschloe commented Aug 31, 2017

pv commented Sep 3, 2017

ENH: Wrap LAPACK's dsytrd #7780

ENH: Wrap LAPACK's dsytrd #7780

Conversation

nschloe commented Aug 22, 2017

ilayn left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nschloe commented Aug 22, 2017

nschloe commented Aug 23, 2017

ev-br commented Aug 23, 2017

nschloe commented Aug 23, 2017

ev-br commented Aug 23, 2017

nschloe commented Aug 24, 2017

ilayn left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ilayn commented Aug 25, 2017

eric-wieser commented Aug 27, 2017

ilayn commented Aug 27, 2017

eric-wieser commented Aug 27, 2017

nschloe commented Aug 27, 2017

nschloe commented Aug 28, 2017

nschloe commented Aug 28, 2017 • edited Loading

nschloe commented Aug 28, 2017

eric-wieser commented Aug 28, 2017

nschloe commented Aug 31, 2017

ilayn commented Aug 31, 2017

nschloe commented Aug 31, 2017

pv commented Sep 3, 2017

nschloe commented Aug 28, 2017 •

edited

Loading