Cookbook: Adding a routine to Slycot

Before starting be sure to have read the manual for f2py. The keywords used by f2py will be marked, e.g intent.

When adding a new routine to Slycot three of the following files will need to be modified:

slycot/__init__.py
slycot/analysis.py
slycot/math.py
slycot/synthesis.py
slycot/transform.py 
slycot/src/analysis.pyf
slycot/src/math.pyf
slycot/src/synthesis.pyf
slycot/src/transform.pyf

For example, in case we were to add routine SB01BD we would need to consider only

slycot/__init__.py
slycot/synthesis.py
slycot/src/synthesis.pyf

Specifically we need to add a signature for f2py in synthesis.pyf, write an high level python wrapper in synthesis.py and finally add the relevant import statements in __init__.py. We will review those steps, as performed in my Linux system, in the following discussion.

Creating the signature

Start by calling f2py to create an initial dummy signature

$cd slycot/src
$f2py -h new.pyf SB01BD.f

this will create the following signature inside new.pyf

subroutine sb01bd(dico,n,m,np,alpha,a,lda,b,ldb,wr,wi,nfp,nap,nup,f,ldf,z,ldz,tol,dwork,ldwork,iwarn,info) ! in SB01BD.f
    character :: dico
    integer :: n
    integer :: m
    integer :: np
    double precision :: alpha
    double precision dimension(lda,*) :: a
    integer optional,check(shape(a,0)==lda),depend(a) :: lda=shape(a,0)
    double precision dimension(ldb,*) :: b
    integer optional,check(shape(b,0)==ldb),depend(b) :: ldb=shape(b,0)
    double precision dimension(*) :: wr
    double precision dimension(*) :: wi
    integer :: nfp
    integer :: nap
    integer :: nup
    double precision dimension(ldf,*) :: f
    integer optional,check(shape(f,0)==ldf),depend(f) :: ldf=shape(f,0)
    double precision dimension(ldz,*) :: z
    integer optional,check(shape(z,0)==ldz),depend(z) :: ldz=shape(z,0)
    double precision :: tol
    double precision dimension(*) :: dwork
    integer :: ldwork
    integer :: iwarn
    integer :: info
end subroutine sb01bd

NOTE: After starting this wiki i added the signature of this routine to the repository. If you want to try to perform this procedure yourself you still can. Simply change the signature in the following way:

subroutine my_sb01bd(dico,n,m,np,alpha,a,lda,b,ldb,wr,wi,nfp,nap,nup,f,ldf,z,ldz,tol,dwork,ldwork,iwarn,info) ! in SB01BD.f
    fortranname sb01bd
    [...cut...]
end subroutine my_sb01bd

this will create an additional wrapper for SB01BD named my_sb01bd.

We can copy and paste it into slycot/src/synthesis.pyf and start editing it. First let’s add some checks to the integers n , m and np

    integer check(n>=0) :: n
    integer check(m>=0) :: m
    integer check(np>=0 && np<=n),depend(n) :: np

here we check for non-negativity, additionally we require np, the number of eigenvalues to be assigned, to be less or equal than the system dimension (&& stands for logical AND in C). Note that since we use the value of n in the line for np we need to make sure that n is initialized before np, this is accomplished by using the keyword depend. Since there is no default value they will be treated as inputs, so the keyword intent or required should not be necessary.

For input alpha, following the routine description we would need to check for positivity whenever dico= ‘D’, something like

double precision check(dico == 'C' || alpha>0),depend(dico):: alpha

where || stands for logical OR in C. Unfortunately this does not work as the first comparison cannot be carried by f2py (or better yet i didn’t find a way to do so). In some other cases this led me to produce more than one f2py wrappers for a single routine. In this case i deem it unnecessary as most probably the routines has an internal check and nothing else seems affected.

The part concerning the system dynamics matrix a could we modified as the following

    double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
    integer intent(hide),depend(a) :: lda=shape(a,0)

here a is both an input and an output thus the intent statement, copy may or may not be used. If it isn’t the input content of a will be overwritten by the routine. The dimension of input arguments need not to be specified in general but it is useful doing so, in fact the wrapper will check if the actual dimension of the matrix are greater or equal than the specified dimensions and throw an eventual error. The routine description will provide the needed minimum dimensions. Note that lda is completely specified by a thus it can be safely hidden.

Similarly we can modify the lines concerning b , wr and wi

    double precision dimension(n,m),depend(n,m) :: b
    integer intent(hide),depend(b) :: ldb=shape(b,0)
    double precision intent(in,out,copy),dimension(np),depend(np) :: wr
    double precision intent(in,out,copy),dimension(np),depend(np) :: wi

The integer outputs can be handled easily

    integer intent(out) :: nfp
    integer intent(out) :: nap
    integer intent(out) :: nup

while the matrix outputs are modified similarly as they were inputs

    double precision intent(out),dimension(m,n) :: f
    integer intent(hide),depend(f) :: ldf=shape(f,0)
    double precision intent(out),dimension(n,n) :: z
    integer intent(hide),depend(z) :: ldz=shape(z,0)

note that the specification of the dimensions for a matrix output is required as the wrapper itself will create them.

The remaining variables are common among SLICOT routines and should be handled consistently. tol is made optional with default value 0.0, this will induce the routine to use its own default value for the tolerance.

double precision :: tol = 0

dwork and ldwork are used as floating point cache, ldwork is made optional with default vale to the minimum required (or some close vale if the computation is too verbose or impossible). In this case the minimum value specified by the routine is

max(1,5*m,5*n,2*n+4*m)

note that max in the signature accept only two values. Two alternatives may be either

    integer :: ldwork = max(5*m,max(5*n,2*n+4*m))+1

or

    integer :: ldwork = max(5*n,2*n+5*m)+1

depending on how tight we want the value.

The line about dwork should be standard among different signatures

     double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork

this tells f2py to use an internal cache array common to all the routines that is expanded as necessary.

Finally the error/warning codes, simply

    integer intent(out) :: iwarn
    integer intent(out) :: info

This lead to the signature

subroutine sb01bd(dico,n,m,np,alpha,a,lda,b,ldb,wr,wi,nfp,nap,nup,f,ldf,z,ldz,tol,dwork,ldwork,iwarn,info) ! in SB01BD.f
    character :: dico
    integer required,check(n>=0) :: n
    integer required,check(m>=0) :: m
    integer required,check(np>=0 && np<=n),depend(n) :: np
    double precision :: alpha
    double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
    integer intent(hide),depend(a) :: lda=shape(a,0)
    double precision dimension(n,m),depend(n,m) :: b
    integer intent(hide),depend(b) :: ldb=shape(b,0)
    double precision intent(in,out,copy),dimension(np),depend(np) :: wr
    double precision intent(in,out,copy),dimension(np),depend(np) :: wi
    integer intent(out) :: nfp
    integer intent(out) :: nap
    integer intent(out) :: nup
    double precision intent(out),dimension(m,n) :: f
    integer intent(hide),depend(f) :: ldf=shape(f,0)
    double precision intent(out),dimension(n,n) :: z
    integer intent(hide),depend(z) :: ldz=shape(z,0)
    double precision :: tol = 0
    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
    integer :: ldwork
    integer intent(out) :: iwarn
    integer intent(out) :: info
end subroutine sb01bd

We can now try to compile and try out the low level wrapper, first touch the file slycot/src/_wrapper.pyf by doing

$touch _wrapper.pyf

this will change the last-modified time and induce f2py to recompile the wrappers. Then the usual setup command from the project root directory

#python setup.py install

Then in IPython we can check if the signature that we wrote make sense

In [1]: import slycot; print slycot._wrapper.sb01bd.__doc__ 
sb01bd - Function signature:
  a,wr,wi,nfp,nap,nup,f,z,iwarn,info = sb01bd(dico,n,m,np,alpha,a,b,wr,wi,[tol,ldwork,overwrite_a,overwrite_wr,overwrite_wi])
Required arguments:
  dico : input string(len=1)
  n : input int
  m : input int
  np : input int
  alpha : input float
  a : input rank-2 array('d') with bounds (n,n)
  b : input rank-2 array('d') with bounds (n,m)
  wr : input rank-1 array('d') with bounds (np)
  wi : input rank-1 array('d') with bounds (np)
Optional arguments:
  overwrite_a := 0 input int
  overwrite_wr := 0 input int
  overwrite_wi := 0 input int
  tol := 0 input float
  ldwork := max(5*n,2*n+5*m)+1 input int
Return objects:
  a : rank-2 array('d') with bounds (n,n)
  wr : rank-1 array('d') with bounds (np)
  wi : rank-1 array('d') with bounds (np)
  nfp : int
  nap : int
  nup : int
  f : rank-2 array('d') with bounds (m,n)
  z : rank-2 array('d') with bounds (n,n)
  iwarn : int
  info : int

Note the optional boolean inputs overwrite_a , overwrite_wr and overwrite_wi added by f2py since we used the intent copy. If 1 no copy is made and both the input and output references of each of a , wr and wi will point to the same memory area, i.e the input values will be overwritten. Default is 0.