Skip to content

Update Slycot to use LAPACK 3.6.x + #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 22 commits into from
Jan 13, 2018
Merged

Update Slycot to use LAPACK 3.6.x + #13

merged 22 commits into from
Jan 13, 2018

Conversation

repagh
Copy link
Member

@repagh repagh commented Jun 10, 2017

Fixed Slycot for newer LAPACK, replaced deprecated function use

repagh and others added 5 commits June 10, 2017 17:50
- add SLCT_DLATZM.f, to replace removed LAPACK function
- adapt various routines to use
  * replacement SLCT_DLATZM
  * DGEGS by DGGES (and work vector size changed)
  * ZLATZM by ZUNMRZ
return locations where result value was not set
Copy link
Collaborator

@roryyorke roryyorke left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks very much for working on this.

This is a superficial review: I haven't studied the Fortran changes in detail, and it'll take me some time to understand them.

Ideally we'd add tests to check the correctness of these. Slycot testing is so sparse I doubt all of these changes (or any!) are covered.

It would also be good to add a test that exercises this change properly; e.g., if I understamd @jgoppert , he hit this problem when linking against lapack 3.6 provided by conda-forge.

Here and there is trailing whitespace, e.g., AB13AD.f line 200 (end of line with "Jun. 2017 ..") (I only noticed because Emacs highlights trailing whitespace in the diff).

@@ -7,7 +7,7 @@ build:

script:
- cd $RECIPE_DIR/..
- $PYTHON setup.py install
- LDFLAGS="-shared" FFLAGS=-fPIC $PYTHON setup.py install # [not win]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you know why this is required? "-shared" and "-fPIC" are standard enough that I would think the build tools we're using (distutils, setuptools, conad) would set them.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't know why, but it would not "conda build" build without, I tried both on OSX and Linux. Tomorrow I can fire up a virtualbox with windows, and see what conda build does there

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Once this works, I will try to get slycot + control in conda-forge

There might be something in the build tools that is not set right, but I don't know where to look for that.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've removed this and the setup.py mods from your branch here; passing Travis CI results here.

C
AB13AD = ZERO
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this related to the deprecated LAPACK functions? If not, could we maybe do these in a separate PR?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no, this was just to get the warnings out.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't do that in a separate branch though. Would be some work to split it off.

@@ -199,8 +199,8 @@ SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD,
INTEGER ILAENV
EXTERNAL ILAENV
C .. External Subroutines ..
EXTERNAL DLAPMT, DLARFG, DLASET, DLATZM, DORMQR, DORMRQ,
$ MB03OY, MB03PY, XERBLA
EXTERNAL DLAPMT, DLARFG, DLASET, SLCT_DLATZM, DORMQR,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't have a particularly strong view on this, but why not just call it DLATZM? In my admittedly limited and rusty understanding of linking, this should work whether or not the LAPACK linked against provides DLATZM.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also possible. I came across this in the LAPACK license:

Like all software, it is copyrighted. It is not trademarked, but we do ask the following:

If you modify the source for these routines we ask that you change the name of the routine and comment the changes made to the original.

We will gladly answer any questions regarding the software. If a modification is done, however, it is the responsibility of the person who modified the routine to provide support.

This makes it clearer we are using our own copy. The six-letter limit for fortran variable names is only relevant for FORTRAN 4 or 5, I think it was already dropped for FORTRAN 77

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, fair enough.

FWIW, my go-to reference for Fortran 77 says "One of the most irksome restrictions of Fortran77 is that symbolic names cannot be more than six characters long."

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK then, our machine had an extension to its FORTRAN 77 then. I actually sold my FORTRAN book to some Indonesian visitors years ago ....

slycot/setup.py Outdated
abiflags = sys.abiflags
except AttributeError:
abiflags = ''
liblist = ['lapack', 'python'+pyver+abiflags]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

similarly to my comment on meta.yaml: why is this needed now?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Linking on Mac OSX failed; when linking a shared library on mac OSX (and on some strict Linux platforms as well, SUSE is one of the stricter ones), the presence of all symbols that the library in turn requires is checked. Only lapack was in the list, but for the python binding also libpython is required. The alternative is relaxing the linking with flags, but this is much better, you at least also know that all symbols are present, and not figure that out when you try to load the module.

CALL ZLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1,
$ DCONJG( TC ), ABCD(IROW,I1+1),
$ ABCD(IROW+1,I1+1), LDABCD, ZWORK )
C CALL zlatzm( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My opinion is it would be better to remove code, rather than comment it out. Having said that, I see you have some questions, so I guess this is work-in-progress.

@@ -550,7 +550,7 @@ subroutine sg03ad(dico,job,fact,trans,uplo,n,a,lda,e,lde,q,ldq,z,ldz,x,ldx,scale
double precision intent(out), dimension(n), depend(n) :: beta
integer intent(hide,cache),dimension(n*n),depend(n) :: iwork
double precision intent(hide,cache),dimension(ldwork) :: dwork
integer optional,check(ldwork>=max(1,n)),depend(n) :: ldwork=max(1,max(2*n*n,4*n))
integer optional,check(ldwork>=max(1,max(2*n*n,8*n+16))),depend(n) :: ldwork=max(1,max(2*n*n,8*n+16))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just to check: is this a bug, or due to the change from DGEGS to DGGES?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

due to the change, I hit one of these (don't remember which) on running the tests, and then I updated & checked the work arrays required on all the DGGES using code. 8*n+16 is the new number of required work doubles.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But now I see it, the check formula was a bug too. Would hit you if you call the wrapper with the work array provided and too small (otherwise it is created by the wrapper with the size calculated according to the 2nd formula).

@repagh
Copy link
Member Author

repagh commented Jun 19, 2017

I removed the trailing whitespace, cleaned up comment, and added a test case for one of the changed routines (sg03ad).

@repagh
Copy link
Member Author

repagh commented Jun 19, 2017

added another test case, for sg02ad

Can someone pull it in? After that I want to try getting slycot and control into conda-forge.

@roryyorke
Copy link
Collaborator

I'm still looking at the Fortran changes you've made; I'm afraid my Fortran was never that good, and I also find the LAPACK API difficult to work with.

Three things puzzle me: first, it looks to me like in SLICOT DLATZM is typically used after DLARFG. I tried to use a simple 3x3 example to explore this, and where I could give DLATZM a 2-element argument representing the Householder transformation vector (this is V in the linked docs), DORMRZ always needed 3-element vector (A in the linked docs), though it was important to set K to 2. This seems odd, but it looks like DORMRZ was designed for use with DTZRZF, which sets up a vector (or, in general, a matrix) of the type DORMRZ expects.

The second thing is that you said SLICOT seems to use DLATZM in a non-standard way, but here's a call to ZLATZM (from AB8NX.f, line 316, a618aeb)

               CALL ZLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1,
     $                      TC )
               CALL ZLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1,
     $                      DCONJG( TC ), ABCD(IROW,I1+1),
     $                      ABCD(IROW+1,I1+1), LDABCD, ZWORK )

and here's a call to DLATZM (AB08NX.f, line 310, also a618aeb):

               CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T )
               CALL DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T,
     $                      ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD,
     $                      DWORK )

These look fairly analogous, so it's odd that you could use ZUNMRZ in the first and not DORMRZ in the second.

The third thing that I haven't looked at properly is that ZUNMRZ (and DORMRZ) don't obviously need larger working areas than the deprecated functions, so I don't understand why that had to change.

In terms of conda-forge: wouldn't it be possible to start off pointing the conda-forge repository at your branch, especially during an initial experimental stage? Presumably that could be changed later.

@repagh
Copy link
Member Author

repagh commented Jun 20, 2017

Looking at it again, it might be possible at several places, but not all.
DLATZM and ZLATZM return the result matrix C in two different variables, C1 for the first row or column and C2, for the remaining rows or columns. (ZUNMRZ and DORMRZ have just one in/out variable, and use/return the complete matrix, without potentially splitting this up.

Substitution by ZUNMRZ / DORMRZ without copying then works only when:

  • in case side = 'L', C1 gives the first row, C2 remaining, and [C1;C2] together are a single contiguous matrix C.
  • in case side ='R', C1 gives first column, C2 remaining columns, and [C1 C2] together is a single contiguous matrix.

In the examples above the counting is explicit (IROW and IROW+1) so it is easily verified. In SB06ND the code uses JKCUR and JCUR as index, and treating input/output as one matrix would only be possible if JCUR is JKCUR+1

               DO 80 KK = KMAX, KMAX - KSTEP + 1, -1
               KCUR = KSTAIR(KK)
               JCUR = JKCUR - KCUR
C
               DO 60 J = 1, KCUR
                  CALL SLCT_DLATZM( 'Left', KCUR+1, N-JCUR+1,
     $                         F(1,JKCUR), 1,
     $                         DWORK(JKCUR), A(JKCUR,JCUR),
     $                         A(JCUR,JCUR), LDA, DWORK(N+1) )
                  JCUR = JCUR - 1
                  JKCUR = JKCUR - 1
   60          CONTINUE
C
   80       CONTINUE

I am guessing, from the code, that that is not the case, KSTAIR would have to have -1 in all cases then. I think this is a piece of code that exploits the properties of DLATZM in a way that is not replaceable by DORMRZ

@roryyorke
Copy link
Collaborator

Thanks for explaining. That is unfortunate! For SIDE='L', one could perhaps use DLASWP to swap rows before and after the DLATZM call (if JCUR>=2), but that doesn't help for SIDE='R', and also feels quite hacky.

What do you think of my first point (ZUNMRZ wanting a matrix like one returned by ZTZRZF) ? I'm not entirely sure it's related, but I've compared the SLICOT example for AB08NZ, which uses AB8NXZ, with "master" (current Slycot master) and your ("rvp") versions, and the results are different. If I run "make test" on the code in the attached archive, I see:

diff -u master.res rvp.res 
--- master.res	2017-06-21 21:56:45.590914206 +0200
+++ rvp.res	2017-06-21 21:56:45.594914206 +0200
@@ -22,15 +22,7 @@
    -4.0000  +0.0000i 
 
 
- The number of finite invariant zeros =   2
-
- The finite invariant zeros are 
-
- real  part     imag  part 
-    2.0000  +0.0000i 
-   -1.0000  +0.0000i 
-
- which correspond to the generalized eigenvalues of (lambda*BF - AF).
+ The number of finite invariant zeros =   0
 
 
  The number of infinite zeros =   2
@@ -43,4 +35,4 @@
  The number of left Kronecker indices =   1
 
  The left Kronecker (row) indices of (A,B,C,D) are 
-  2
+  4

This is on an Ubuntu 14.04 x64 system. The master results are identical to the reference results.

I got the example from the Debian source package of SLICOT.

If I'm right (and I have no idea if I am), maybe we should instead go with a SLCT_ZLATZM?

@repagh
Copy link
Member Author

repagh commented Jun 24, 2017

I reverted to SLCT_ZLATZM, does that fix the test?

Checked. It does seem to fix it.

Thanks for looking up the Debian slicot. It seems a bit newer than the slicot we are using, there were some routines added. The set of examples/tests is also quite good. Our SB04OD.f has a segfault which is fixed in the debian version. Should fix that too, but in another request?

@roryyorke
Copy link
Collaborator

Are you going to use the new conda-forge package, or do you still require this change for your course?

I think it would be good to sort out the {D,Z}LATZM deprecation, but if there's no hurry I'd like to more carefully review the build changes, etc., proposed here, and perhaps add some tests along the lines of the example for AB08NZ.

It would be good to update to the version of SLICOT Debian has - perhaps we should do that (in a separate PR) before trying to fix the deprecation problems.

@roryyorke
Copy link
Collaborator

A few questions; apologies if they have obvious answers

  1. Should we remove our own conda recipes? They don't seem useful anymore; we could just use the usual "python setup.py build" to build and test on Travis-CI, and for development on our own systems.
  2. When doing testing on Travis-CI, should python-control install Slycot from conda-forge?
  3. Is there any reason to prefer OpenBLAS to MKL for Slycot? MKL seems to be the standard library installed with numpy when using conda, so it might be a good idea to save users from downloading OpenBLAS if possible.
  4. Is there any reason not to make python-control a "noarch" package?

@repagh
Copy link
Member Author

repagh commented Jul 19, 2017

I am currently also looking at the new conda-forge package for my control course. Linux and Mac OSX versions work. However, most students still use Windows, and there is no Windows version. I have tried both approaches, building windows packages with the conda-forge recipe, however that requires openblas, and that is not available on windows. I also tried to build openblas, but I for now hit a dead end there.

My own conda packages (repa / slycot) work with lapack and blas.

I am also trying to build the new slycot (my own repository for now), because blas + lapack are available for Windows. But I am still struggling & stuck at the linking process, although it seems it is just a question of finding the right libraries to link. I am starting to learn more than I ever cared to know about compiling python extensions on Windows.

I don't have answers to your questions.

@repagh
Copy link
Member Author

repagh commented Oct 9, 2017

I finally finally managed to create a build for Windows (64 bit) too. It is on my conda channel now.
To pull that off, I had to use mingw32 build. I tried all sorts of gfortran + msvc combinations, which involved hunting down libraries from both the mingw and msvc domains, but all were fruitless.

  • The bld.bat now selects a mingw32 build
  • This uses the lapack build for windows from conda-forge (with conda-forge as a secondary channel)
  • Numpy needs to be 1.13.x . The new python interpreters from anaconda are built with msvc 2015, (numbered 1900 / 140 ????), and older numpy distutils have no configuration values for this version. The mingw32 build checks the msvc version, and fails for older numpy. For some unknown reason, I had to actually force numpy to the latest versions in meta.yaml.
  • However, numpy 1.13.x is not available for py33 or py34 for anaconda, and these two therefore now fail on Travis.

The conda recipe, and the build method, are different from the ones in conda-forge. This is apparently the only slycot version building for new anaconda in Windows; until Windows building works with the conda-forge versions, I would keep the recipes in.

I don't specify a blas version in the meta.yaml, so the build will pull in what lapack and numpy need.

@ilayn
Copy link

ilayn commented Oct 10, 2017

@repagh As a wonderful achievement of the wizards, SciPy wheels finally would be available for windows too in version 1.0 which would make it possible for everyone to use pip to install.

This is done with a magical combination of mingw and msvc. The .yml configuration for Appveyor and also the MacPython version might guide you to automate the build for windows and if it works probably will solve the Slycot torture for windows users in general via pip or conda installation.

@murrayrm
Copy link
Member

murrayrm commented Dec 24, 2017

I've done some testing of this PR against python-control and the new version appears to work against lapack3.5 (which is the default for Travis CI, which uses Ubuntu Testy). Details available on my local build.

I want to do some testing to make sure this version works on MacOS and against the proper versions of lapack (3.6+), but if all of that works then I will merge it in to master.

@murrayrm
Copy link
Member

I've done some testing on my local machine (MacOS 10.12.6) and I was able to get everything to compile and (basically) run correctly. For posterity, I am using the following commands to build and install slycot:

git checkout -b repagh-master master
git pull https://github.com/repagh/Slycot.git master
conda create -n test_slycot -c conda-forge python=3.6 gcc lapack matplotlib scipy
source activate test_slycot
python setup.py install

This gives the following versions of key packages (using conda list):

gcc                       4.8.5                         8  
lapack                    3.6.1                         1    conda-forge
python                    3.6.4                         0    conda-forge
scipy                     1.0.0           py36_blas_openblas_201  [blas_openblas]  conda-forge

(note that scipy turns out not to be the default one, but comes from conda-forge; not sure if this matters.)

I then go over to python-control and run python setup.py test. This gives the following error:

FAIL: testMinrealBrute (control.tests.minreal_test.TestMinreal)
Traceback (most recent call last):
  File "/Users/murray/Dropbox/macosx/src/python-control/murrayrm/control/tests/minreal_test.py", line 62, in testMinrealBrute
    raise e
  File "/Users/murray/Dropbox/macosx/src/python-control/murrayrm/control/tests/minreal_test.py", line 57, in testMinrealBrute
    ht1.den[0][0], ht2.den[0][0])
  File "/Users/murray/Dropbox/macosx/src/python-control/murrayrm/control/tests/minreal_test.py", line 35, in assert_numden_almost_equal
    np.testing.assert_array_almost_equal(n1, n2)
  File "/Users/murray/anaconda/envs/test_slycot/lib/python3.6/site-packages/numpy/testing/utils.py", line 962, in assert_array_almost_equal
    precision=decimal)
  File "/Users/murray/anaconda/envs/test_slycot/lib/python3.6/site-packages/numpy/testing/utils.py", line 715, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Arrays are not almost equal to 6 decimals

(shapes (3,), (4,) mismatch)
 x: array([ -1.33214 , -10.368378, -13.703545])
 y: array([ -1.33214 , -13.297947, -36.505109, -30.136079])

I'm not sure what this error is about, but I get the exact same error when I use the current slycot master => probably not something that is wrong in this PR. All other tests are OK:

Ran 247 tests in 9.653s

FAILED (failures=1, skipped=2)
Test failed: <unittest.runner.TextTestResult run=247 errors=0 failures=1>
error: Test failed: <unittest.runner.TextTestResult run=247 errors=0 failures=1>

As noted above, the python-control tests also seem to work correctly in Travis CI => something different about the Mac version, but doesn't seem to be limited to the new version in this PR.

Based on this testing, I think it is OK to merge this into master. In particular, I note the following:

  • This PR works correctly using the old version of lapack (3.5) under Ubuntu/Testy (via Travis CI)
  • This PR works correctly using the new version of lapack (3.6) under MacOS/Sierra (modulo the errors above)
  • This PR works correctly in various environments tested above by @repagh.

@repagh and @roryyorke: since the two of you have looked at this in more detail, it would be good to get your OK on the testing above as sufficient for merging.

@murrayrm
Copy link
Member

@roryyorke Regarding your questions:

Should we remove our own conda recipes? They don't seem useful anymore; we could just use the usual "python setup.py build" to build and test on Travis-CI, and for development on our own systems.

Not sure whether the conda recipes are being used by anyone, but I note that in python control PR 169 I have removed the use of conda-build from Travis CI => we no longer need the conda recipe for Travis CI.

When doing testing on Travis-CI, should python-control install Slycot from conda-forge?

I have removed this functionality in python control PR 169. The new .travis.yml file builds slycot from scratch => no longer reliant on conda-forge.

Is there any reason to prefer OpenBLAS to MKL for Slycot? MKL seems to be the standard library installed with numpy when using conda, so it might be a good idea to save users from downloading OpenBLAS if possible.

No strong preference here, but things have worked OK for me using lapack from conda-forge (see comments above).

Is there any reason not to make python-control a "noarch" package?

This discussion probably belongs in python-control?

@roryyorke
Copy link
Collaborator

I've had another quick look over this, and it looks fine to me. Having "own" copies of the deprecated LAPACK functions is the right solution given the limited resources we have.

@repagh has addressed my review comments.

I have not reviewed the Windows build changes (9ce23a8 and later), and can't comment on those.

Thanks for the response to my other questions. I haven't looked at python-control/python-control#169 yet. Two things occur to me:

  1. This will reduce commit-result turnaround time. Possibly not that big a deal.
  2. Will python-control always then build against HEAD of Slycot?

Thanks for looking at this. I'm afraid I don't have much free time at the moment, but I'll try to look at the other PRs you've made.

@murrayrm murrayrm merged commit a41ac6e into python-control:master Jan 13, 2018
bnavigator pushed a commit to bnavigator/SLICOT-Reference that referenced this pull request Feb 7, 2021
Ported from python-control/Slycot#13: commits 2d9a253 and 8cf5e49

Includes c1cb581 add an explicit return of ZERO in functions that have alternative
bnavigator added a commit to bnavigator/SLICOT-Reference that referenced this pull request Feb 7, 2021
bnavigator added a commit to bnavigator/SLICOT-Reference that referenced this pull request Feb 7, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants