I abandoned the ctypes route to exposing dtgsen and am trying the f2py route.  I'm hoping someone might be able to answer my questions.  The f2py documentation is pretty good, and I am beginning to feel like I have some idea of what I am doing, but I am running into problems with a few f2py directives.  For the record, I am only interested in the function for the case when IJOB=0. 

I have produced a file ordqz.so that I can import and look at the docstring, but when I try to use the function I get  

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

Which I assume comes from the LWORK argument as stated in the docs.

I have defined M, PL, PR, WORK, LWORK, IWORK, and LIWORK as follows in ordqz.pyf based on the docs <http://linux.die.net/man/l/dtgsen> (in-line below).  These should be defined based on the value of IJOB.  Is there some other way to go about this, or can someone just point me to a relevant function that's already in the numpy/scipy source as a reference for a similar situation?

!xxx - made up dimensions for pl and pr, and dependency of m
           integer depend(pl),intent(hide):: m = shape(pl,0)
           double precision intent(out),dimension(n):: pl
           double precision intent(out),dimension(n):: pr
           double precision dimension(2),intent(out):: dif
           double precision intent(hide,cache),dimension(lwork),depend(lwork):: work
!Made up a value for lwork (its lower bound) 4*n+16
           integer intent(hide),depend(n):: lwork=4*n+16
           integer intent(hide,cache),dimension(liwork),depend(liwork):: iwork
!Made up a value for liwork (its lower bound) n+6
           integer intent(hide):: liwork=n+6


Relevant part of docs
-------
M (output) INTEGER
   The dimension of the specified pair of left and right eigen- spaces (deflating subspaces). 0 <= M <= N.

   PL, PR (output) DOUBLE PRECISION If IJOB = 1, 4 or 5, PL, PR are lower bounds on the reciprocal of the norm of "projections" onto left and right eigenspaces with respect to the selected cluster. 0 < PL, PR <= 1. If M = 0 or M = N, PL = PR = 1. If IJOB = 0, 2 or 3, PL and PR are not referenced.
DIF (output) DOUBLE PRECISION array, dimension (2).
   If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
   If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
   Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based estimates of Difu and Difl. If M = 0 or N, DIF(1:2) = F-norm([A, B]). If IJOB = 0 or 1, DIF is not referenced.
WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
   IF IJOB = 0, WORK is not referenced. Otherwise, on exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK (input) INTEGER
   The dimension of the array WORK. LWORK >= 4*N+16. If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).

   If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.
IWORK (workspace/output) INTEGER array, dimension (LIWORK)
   IF IJOB = 0, IWORK is not referenced. Otherwise, on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
LIWORK (input) INTEGER
   The dimension of the array IWORK. LIWORK >= 1. If IJOB = 1, 2 or 4, LIWORK >= N+6. If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6).

   If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA.

Thanks for any help,

Skipper