Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

ENH: Add QZ decomposition with tests #185

Closed
wants to merge 8 commits into from

4 participants

@jseabold
Collaborator

This adds the generalized Schur decomposition to scipy.linalg. It has been requested on the mailing list a few times. I didn't follow the templating pattern in flapack_user because I have no idea what the syntax for the generic signature in gees__user__routines means.

scipy/linalg/decomp_qz.py
((146 lines not shown))
+ typb = 'F'
+
+ overwrite_a = overwrite_a or (_datacopied(a1,A))
+ overwrite_b = overwrite_b or (_datacopied(b1,B))
+
+ gges, = get_lapack_funcs(('gges',), (a1,b1))
+
+ if lwork is None or lwork == -1:
+ # get optimal work array size
+ result = gges(lambda x: None, a1, b1, lwork=-1)
+ lwork = result[-2][0].real.astype(np.int)
+
+ if sort is None:
+ sort_t = 0
+ sfunction = lambda x : None
+ else: #TODO, how to specify the external user function?
@pv Owner
pv added a note

Is this comment outdated?

@jseabold Collaborator

yes, removed

@jseabold Collaborator

Although, I still don't know how to use the templated external functions as what was already in scipy/linalg/flapack_user.pyf.src. I have no idea what the call signature here means or how to add an optional third parameter as would be needed.

function <prefix>select(<_arg=arg1\,arg2,\0,arg,\2>)

I pinged the dev list, but no one replied. But the hard-coded ones I added work.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/linalg/flapack.pyf.src
((33 lines not shown))
+ integer optional,intent(in),depend(jobvsr,n) :: ldvsr=((jobvsr==1)?n:1)
+ <ftype2> intent(out),depend(lwork),dimension(MAX(lwork,1)) :: work
+ integer optional,intent(in),depend(n),check(lwork>=MAX(1,8*n+16)||lwork==-1):: lwork=8*n+16
+ logical optional,intent(hide),depend(n),dimension(n) :: bwork
+ integer intent(out):: info
+ end subroutine <prefix2>gges
+
+ subroutine <prefix2c>gges(jobvsl,jobvsr,sort_t,<prefix2c>select,n,a,lda,b,ldb,sdim,alpha,beta,vsl,ldvsl,vsr,ldvsr,work,lwork,rwork,bwork,info)
+
+ ! For a pair of N-by-N complex nonsymmetric matrices (A,B) computes
+ ! the generalized eigenvalues, the generalized real Schur form (S,T),
+ ! optionally, the left and/or right matrices of Schur vectors (VSL
+ ! and VSR).
+ ! (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**H )
+
+ callstatement (*f2py_func)((jobvsl?"V":"N"),(jobvsr?"V":"N"),(sort_t?"S":"N"),cb_<prefix2c>select_in_<prefix2c>gges__user__routines,&n,a,&lda,b,&ldb,&sdim,alpha,beta,vsl,&ldvsl,vsr,&ldvsr,work,&lwork,rwork,bwork,&info,1,1)
@pv Owner
pv added a note

What are the two parameters after &info ? They're not there in the *gges manpage...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/linalg/flapack.pyf.src
@@ -15,6 +16,78 @@ interface
include 'flapack_user.pyf.src'
+ subroutine <prefix2>gges(jobvsl,jobvsr,sort_t,<prefix2>select,n,a,lda,b,ldb,sdim,alphar,alphai,beta,vsl,ldvsl,vsr,ldvsr,work,lwork,bwork,info)
+
+ ! For a pair of N-by-N real nonsymmetric matrices (A,B) computes
+ ! the generalized eigenvalues, the generalized real Schur form (S,T),
+ ! optionally, the left and/or right matrices of Schur vectors (VSL
+ ! and VSR).
+ ! (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
+
+ callstatement (*f2py_func)((jobvsl?"V":"N"),(jobvsr?"V":"N"),(sort_t?"S":"N"),cb_<prefix2>select_in_<prefix2>gges__user__routines,&n,a,&lda,b,&ldb,&sdim,alphar,alphai,beta,vsl,&ldvsl,vsr,&ldvsr,work,&lwork,bwork,&info,1,1,1)
@pv Owner
pv added a note

What are the parameters after &info?

@jseabold Collaborator

I have no idea. I followed the wrapper for *gees and it has these extraneous parameters as well AFAICT. I assumed it had something to do with the external function call since in that case select it takes at most two parameters.

@pv Owner
pv added a note

I'm starting to believe that the *gees declaration is in error: the two arguments are just tacked on to the Fortran routine prototype, so the routine gets called with the wrong signature. Apparently, this does no harm...

I'd suggest you just remove the two last arguments from the callstatement and callprotoargument. Things should still work OK.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv
Owner
pv commented

This looks good (didn't test it yet, though).

I'd maybe rename the file to _decomp_qz.py to explicitly marked as private, even if the other ones aren't.

@jseabold
Collaborator

Renamed. For some reason I never received a notification on these comments. I just happened to check in.

scipy/linalg/__init__.py
@@ -134,6 +135,7 @@
from decomp_lu import *
from decomp_cholesky import *
from decomp_qr import *
+from decomp_qz import *
@rgommers Owner

Should be from _decomp_qz import *.

@jseabold Collaborator

Right, thanks. Fixed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
Owner

Works fine on OS X 10.6, Python 2.6

@pv pv commented on the diff
scipy/linalg/flapack_user.pyf.src
@@ -6,3 +6,45 @@ python module gees__user__routines
end function <prefix>select
end interface
end python module gees__user__routines
+
+python module zgges__user__routines
@pv Owner
pv added a note

Here you could use f2py's templating:

The syntax <_arg=arg1\,arg2,\0,arg,\2> means the same as <_arg=arg1\,arg2,arg1\,arg2,arg,arg> -- one of the slashes escapes the comma, and apparently the syntax \0, \1, \2, ... can be used (but IMO shouldn't be because it's confusing) to refer to other values in the expansion.

< ftype> :: <_arg>

expands to

real :: arg1,arg2
double precision :: arg1, arg2
complex :: arg
double complex :: arg

I'd maybe also not name the dummy arguments in such a complicated way --- someone reading this code later on will be confused why there are so many underscores, and whether they are needed by f2py for something magical...

@jseabold Collaborator
jseabold added a note

Sorry forgot to look into this. Thanks for the pointer. This is still very murky to me. I still don't see why arg, \2 expands to arg2, arg, arg and why arg, arg should be assumed to be a complex number.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
Owner

Could you also add .. versionadded:: 0.11.0 to the Notes section of the docstring, and mention this in the release notes?

@albop

On my computer Ubuntu 12.04 64, the test test_qz_double fails.
I get the message: error: (lwork>=MAX(1,8*n+16)||lwork==-1) failed for 6th keyword lwork: dgges:lwork=52
Looking at the code, I see that the parameter lwork computed line 156 is indeed incorrect for double types. This can be reproduced by running :

from scipy.linalg.lapack import get_lapack_funcs
n = 5
#A = random([n,n]).astype(complex64)
#B = random([n,n]).astype(complex64)
A = random([n,n]).astype(float64)
B = random([n,n]).astype(float64)
qqz, = get_lapack_funcs(('gges',), (A, B) )
lwork = qqz(lambda x: None, A,B, lwork=-1)[-2][0].astype(int)
print lwork

I have also checked that a direct fortran call with lwork seems to work ( qqz(lambda x: None, A,B) gives the expected result) so I wonder why we need to call the fortran code two times in case the user hasn't supplied lwork to the python function.

Collaborator

Thanks for the report. Unfortunately, I can't reproduce this. What version of lapack are you using?

from scipy.linalg.lapack import get_lapack_funcs
import numpy as np

n = 5
A = np.random.random([n,n]).astype(np.float64)
B = np.random.random([n,n]).astype(np.float64)
qqz, = get_lapack_funcs(('gges',), (A, B))
lwork = qqz(lambda x: None, A,B, lwork=-1)[-2][0].astype(int)
print lwork

lwork = qqz(lambda x: None, A,B)[-2][0].astype(int)
print lwork

Both give 61, though the former call just calculates this value and returns. I have ATLAS 3.9.52 with LAPACK 3.3.1.

Collaborator

Honestly, I have no idea why there are problems. From my end, everything looks ok in the Python code and the Fortran wrappers. I just tried again on another machine that has Ubuntu 10.04 and ATLAS 3.6.0 and LAPACK 3.2.1 from the Ubuntu respositories (my other ones are compiled from source). Everything is fine there even though it gives a different optimal lwork. It is still consistent and works.

I pinged the scipy-dev list to see if anyone else has any insight.

Collaborator

Ok, I found the problem. Could you confirm that the recent changes fix things for you?

Apparently the errors are gone. Good work !

@pv
Owner
pv commented

Merged in 63291d1
Thanks a lot!

@pv pv closed this
@jnothman jnothman referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
@jseabold jseabold deleted the branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
5 doc/release/0.11.0-notes.rst
@@ -67,6 +67,11 @@ A function for creating Pascal matrices, ``scipy.linalg.pascal``, was added.
``misc.logsumexp`` now takes an optional ``axis`` keyword argument.
+QZ Decomposition
+----------------
+
+It is now possible to calculate the QZ, or Generalized Schur, decomposition using `scipy.linalg.qz`. This function wraps the LAPACK routines sgges, dgges, cgges, and zgges.
+
Deprecated features
===================
View
2  scipy/linalg/__init__.py
@@ -68,6 +68,7 @@
cho_solve_banded - Solve previously factored banded linear system
qr - QR decomposition of a matrix
qr_multiply - QR decomposition and multiplication by Q
+ qz - QZ decomposition of a pair of matrices
schur - Schur decomposition of a matrix
rsf2csf - Real to complex Schur form
hessenberg - Hessenberg form of a matrix
@@ -134,6 +135,7 @@
from decomp_lu import *
from decomp_cholesky import *
from decomp_qr import *
+from _decomp_qz import *
from decomp_svd import *
from decomp_schur import *
from matfuncs import *
View
194 scipy/linalg/_decomp_qz.py
@@ -0,0 +1,194 @@
+import warnings
+
+import numpy as np
+from numpy import asarray_chkfinite, single
+
+from misc import LinAlgError, _datacopied
+from lapack import get_lapack_funcs
+
+__all__ = ['qz']
+
+_double_precision = ['i','l','d']
+
+def _select_function(sort, typ):
+ if typ in ['F','D']:
+ if callable(sort):
+ #assume the user knows what they're doing
+ sfunction = sort
+ elif sort == 'lhp':
+ sfunction = lambda x,y: (np.real(x/y) < 0.0)
+ elif sort == 'rhp':
+ sfunction = lambda x,y: (np.real(x/y) >= 0.0)
+ elif sort == 'iuc':
+ sfunction = lambda x,y: (abs(x/y) <= 1.0)
+ elif sort == 'ouc':
+ sfunction = lambda x,y: (abs(x/y) > 1.0)
+ else:
+ raise ValueError("sort parameter must be None, a callable, or "
+ "one of ('lhp','rhp','iuc','ouc')")
+ elif typ in ['f','d']:
+ if callable(sort):
+ #assume the user knows what they're doing
+ sfunction = sort
+ elif sort == 'lhp':
+ sfunction = lambda x,y,z: (np.real((x+y*1j)/z) < 0.0)
+ elif sort == 'rhp':
+ sfunction = lambda x,y,z: (np.real((x+y*1j)/z) >= 0.0)
+ elif sort == 'iuc':
+ sfunction = lambda x,y,z: (abs((x+y*1j)/z) <= 1.0)
+ elif sort == 'ouc':
+ sfunction = lambda x,y,z: (abs((x+y*1j)/z) > 1.0)
+ else:
+ raise ValueError("sort parameter must be None, a callable, or "
+ "one of ('lhp','rhp','iuc','ouc')")
+ else: # to avoid an error later
+ raise ValueError("dtype %s not understood" % typ)
+ return sfunction
+
+
+def qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False,
+ overwrite_b=False):
+ """
+ QZ decompostion for generalized eigenvalues of a pair of matrices.
+
+ The QZ, or generalized Schur, decomposition for a pair of N x N
+ nonsymmetric matrices (A,B) is
+
+ (A,B) = (Q*AA*Z', Q*BB*Z')
+
+ where AA, BB is in generalized Schur form if BB is upper-triangular
+ with non-negative diagonal and AA is upper-triangular, or for real QZ
+ decomposition (output='real') block upper triangular with 1x1
+ and 2x2 blocks. In this case, the 1x1 blocks correpsond to real
+ generalized eigenvalues and 2x2 blocks are 'standardized' by making
+ the correpsonding elements of BB have the form::
+
+ [ a 0 ]
+ [ 0 b ]
+
+ and the pair of correpsonding 2x2 blocks in AA and BB will have a complex
+ conjugate pair of generalized eigenvalues. If (output='complex') or A
+ and B are complex matrices, Z' denotes the conjugate-transpose of Z.
+ Q and Z are unitary matrices.
+
+ Parameters
+ ----------
+ A : array-like, shape (N,N)
+ 2d array to decompose
+ B : array-like, shape (N,N)
+ 2d array to decompose
+ output : str {'real','complex'}
+ Construct the real or complex QZ decomposition for real matrices.
+ lwork : integer, optional
+ Work array size. If None or -1, it is automatically computed.
+ sort : {None, callable, 'lhp', 'rhp', 'iuc', 'ouc'}
+ Specifies whether the upper eigenvalues should be sorted. A callable
+ may be passed that, given a eigenvalue, returns a boolean denoting
+ whether the eigenvalue should be sorted to the top-left (True). For
+ real matrix pairs, the sort function takes three real arguments
+ (alphar, alphai, beta). The eigenvalue x = (alphar + alphai*1j)/beta.
+ For complex matrix pairs or output='complex', the sort function
+ takes two complex arguments (alpha, beta). The eigenvalue
+ x = (alpha/beta).
+ Alternatively, string parameters may be used:
+ 'lhp' Left-hand plane (x.real < 0.0)
+ 'rhp' Right-hand plane (x.real > 0.0)
+ 'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
+ 'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
+ Defaults to None (no sorting).
+
+ Returns
+ -------
+ AA : array, shape (N,N)
+ Generalized Schur form of A.
+ BB : array, shape (N,N)
+ Generalized Schur form of B.
+ Q : array, shape (N,N)
+ The left Schur vectors.
+ Z : array, shape (N,N)
+ The right Schur vectors.
+ sdim : int
+ If sorting was requested, a fifth return value will contain the
+ number of eigenvalues for which the sort condition was True.
+
+ Notes
+ -----
+ Q is transposed versus the equivalent function in Matlab.
+
+ .. versionadded:: 0.11.0
+ """
+ if not output in ['real','complex','r','c']:
+ raise ValueError("argument must be 'real', or 'complex'")
+
+ a1 = asarray_chkfinite(A)
+ b1 = asarray_chkfinite(B)
+
+ a_m, a_n = a1.shape
+ b_m, b_n = b1.shape
+ try:
+ assert a_m == a_n == b_m == b_n
+ except AssertionError:
+ raise ValueError("Array dimensions must be square and agree")
+
+ typa = a1.dtype.char
+ if output in ['complex', 'c'] and typa not in ['F','D']:
+ if typa in _double_precision:
+ a1 = a1.astype('D')
+ typa = 'D'
+ else:
+ a1 = a1.astype('F')
+ typa = 'F'
+ typb = b1.dtype.char
+ if output in ['complex', 'c'] and typb not in ['F','D']:
+ if typb in _double_precision:
+ b1 = b1.astype('D')
+ typb = 'D'
+ else:
+ b1 = b1.astype('F')
+ typb = 'F'
+
+ overwrite_a = overwrite_a or (_datacopied(a1,A))
+ overwrite_b = overwrite_b or (_datacopied(b1,B))
+
+ gges, = get_lapack_funcs(('gges',), (a1,b1))
+
+ if lwork is None or lwork == -1:
+ # get optimal work array size
+ result = gges(lambda x: None, a1, b1, lwork=-1)
+ lwork = result[-2][0].real.astype(np.int)
+
+ if sort is None:
+ sort_t = 0
+ sfunction = lambda x : None
+ else:
+ sort_t = 1
+ sfunction = _select_function(sort, typa)
+
+ result = gges(sfunction, a1, b1, lwork=lwork, overwrite_a=overwrite_a,
+ overwrite_b=overwrite_b, sort_t=sort_t)
+
+ info = result[-1]
+ if info < 0:
+ raise ValueError("Illegal value in argument %d of gges" % -info)
+ elif info > 0 and info <= a_n:
+ warnings.warn("The QZ iteration failed. (a,b) are not in Schur "
+ "form, but ALPHAR(j), ALPHAI(j), and BETA(j) should be correct"
+ "for J=%d,...,N" % info-1, UserWarning)
+ elif info == a_n+1:
+ raise LinAlgError("Something other than QZ iteration failed")
+ elif info == a_n+2:
+ raise LinAlgError("After reordering, roundoff changed values of some"
+ "complex eigenvalues so that leading eigenvalues in the"
+ "Generalized Schur form no longer satisfy sort=True."
+ "This could also be caused due to scaling.")
+ elif info == a_n+3:
+ raise LinAlgError("Reordering failed in <s,d,c,z>tgsen")
+
+ # output for real
+ #AA, BB, sdim, alphar, alphai, beta, vsl, vsr, work, info
+ # output for complex
+ #AA, BB, sdim, alphai, beta, vsl, vsr, work, info
+ if sort_t == 0:
+ return result[0], result[1], result[-4], result[-3]
+ else:
+ return result[0], result[1], result[-4], result[-3], result[2]
View
73 scipy/linalg/flapack.pyf.src
@@ -6,6 +6,7 @@
! $Revision$ $Date$
!
! Additions by Travis Oliphant, Tiziano Zito, Collin RM Stocks, Fabian Pedregosa
+! Skipper Seabold
!
! <prefix2=s,d> <ctype2=float,double> <ftype2=real,double precision>
! <prefix2c=c,z> <ftype2c=complex,double complex> <ctype2c=complex_float,complex_double>
@@ -15,6 +16,78 @@ interface
include 'flapack_user.pyf.src'
+ subroutine <prefix2>gges(jobvsl,jobvsr,sort_t,<prefix2>select,n,a,lda,b,ldb,sdim,alphar,alphai,beta,vsl,ldvsl,vsr,ldvsr,work,lwork,bwork,info)
+
+ ! For a pair of N-by-N real nonsymmetric matrices (A,B) computes
+ ! the generalized eigenvalues, the generalized real Schur form (S,T),
+ ! optionally, the left and/or right matrices of Schur vectors (VSL
+ ! and VSR).
+ ! (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
+
+ callstatement (*f2py_func)((jobvsl?"V":"N"),(jobvsr?"V":"N"),(sort_t?"S":"N"),cb_<prefix2>select_in_<prefix2>gges__user__routines,&n,a,&lda,b,&ldb,&sdim,alphar,alphai,beta,vsl,&ldvsl,vsr,&ldvsr,work,&lwork,bwork,&info)
+ callprotoargument char*,char*,char*,int(*)(<ctype2>*,<ctype2>*,<ctype2>*),int*,<ctype2>*,int*,<ctype2>*,int*,int*,<ctype2>*,<ctype2>*,<ctype2>*,<ctype2>*,int*,<ctype2>*,int*,<ctype2>*,int*,int*,int*
+
+ use <prefix2>gges__user__routines
+
+ integer optional,intent(in),check(jobvsl==0||jobvsl==1) :: jobvsl=1
+ integer optional,intent(in),check(jobvsr==0||jobvsr==1) :: jobvsr=1
+ integer optional,intent(in),check(sort_t==0||sort_t==1) :: sort_t=0
+ external <prefix2>select
+ integer intent(hide),depend(a) :: n=shape(a,1)
+ <ftype2> intent(in,out,copy),dimension(lda,*) :: a
+ integer intent(hide),depend(a) :: lda=shape(a,0)
+ <ftype2> intent(in,out,copy),dimension(ldb,*) :: b
+ integer intent(hide),depend(b) :: ldb=shape(b,0)
+ integer intent(out) :: sdim=0
+ <ftype2> intent(out),dimension(n) :: alphar
+ <ftype2> intent(out),dimension(n) :: alphai
+ <ftype2> intent(out),dimension(n) :: beta
+ <ftype2> intent(out),depend(ldvsl,n),dimension(ldvsl,n) :: vsl
+ integer optional,intent(in),depend(jobvsl,n) :: ldvsl=((jobvsl==1)?n:1)
+ <ftype2> intent(out),depend(ldvsr,n),dimension(ldvsr,n) :: vsr
+ integer optional,intent(in),depend(jobvsr,n) :: ldvsr=((jobvsr==1)?n:1)
+ <ftype2> intent(out),depend(lwork),dimension(MAX(lwork,1)) :: work
+ integer optional,intent(in),depend(n),check(lwork>=MAX(1,MAX(8*n, 6*n+16))||lwork==-1):: lwork=8*n+16
+ logical optional,intent(hide),depend(n),dimension(n) :: bwork
+ integer intent(out):: info
+ end subroutine <prefix2>gges
+
+ subroutine <prefix2c>gges(jobvsl,jobvsr,sort_t,<prefix2c>select,n,a,lda,b,ldb,sdim,alpha,beta,vsl,ldvsl,vsr,ldvsr,work,lwork,rwork,bwork,info)
+
+ ! For a pair of N-by-N complex nonsymmetric matrices (A,B) computes
+ ! the generalized eigenvalues, the generalized real Schur form (S,T),
+ ! optionally, the left and/or right matrices of Schur vectors (VSL
+ ! and VSR).
+ ! (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**H )
+
+ callstatement (*f2py_func)((jobvsl?"V":"N"),(jobvsr?"V":"N"),(sort_t?"S":"N"),cb_<prefix2c>select_in_<prefix2c>gges__user__routines,&n,a,&lda,b,&ldb,&sdim,alpha,beta,vsl,&ldvsl,vsr,&ldvsr,work,&lwork,rwork,bwork,&info)
+ callprotoargument char*,char*,char*,int(*)(<ctype2c>*,<ctype2c>*),int*,<ctype2c>*,int*,<ctype2c>*,int*,int*,<ctype2c>*,<ctype2c>*,<ctype2c>*,int*,<ctype2c>*,int*,<ctype2c>*,int*,<ctype2>*,int*,int*
+
+ use <prefix2c>gges__user__routines
+
+ integer optional,intent(in),check(jobvsl==0||jobvsl==1) :: jobvsl=1
+ integer optional,intent(in),check(jobvsr==0||jobvsr==1) :: jobvsr=1
+ integer optional,intent(in),check(sort_t==0||sort_t==1) :: sort_t=0
+ external <prefix2c>select
+ integer intent(hide),depend(a) :: n=shape(a,1)
+ <ftype2c> intent(in,out,copy),dimension(lda,*) :: a
+ integer intent(hide),depend(a) :: lda=shape(a,0)
+ <ftype2c> intent(in,out,copy),dimension(ldb,*) :: b
+ integer intent(hide),depend(b) :: ldb=shape(b,0)
+ integer intent(out) :: sdim=0
+ <ftype2c> intent(out),dimension(n) :: alpha
+ <ftype2c> intent(out),dimension(n) :: beta
+ <ftype2c> intent(out),depend(ldvsl,n),dimension(ldvsl,n) :: vsl
+ integer optional,intent(in),depend(jobvsl,n) :: ldvsl=((jobvsl==1)?n:1)
+ <ftype2c> intent(out),depend(ldvsr,n),dimension(ldvsr,n) :: vsr
+ integer optional,intent(in),depend(jobvsr,n) :: ldvsr=((jobvsr==1)?n:1)
+ <ftype2c> intent(out),depend(lwork),dimension(MAX(lwork,1)) :: work
+ integer optional,intent(in),depend(n),check(lwork>=MAX(1,2*n)||lwork==-1):: lwork=2*n
+ <ftype2> intent(hide),dimension(8*n) :: rwork
+ logical optional,intent(hide),depend(n),dimension(n) :: bwork
+ integer intent(out):: info
+ end subroutine <prefix2c>gges
+
subroutine <prefix2>pbtrf(lower,n,kd,ab,ldab,info)
! Compute Cholesky decomposition of banded symmetric positive definite
View
42 scipy/linalg/flapack_user.pyf.src
@@ -6,3 +6,45 @@ python module gees__user__routines
end function <prefix>select
end interface
end python module gees__user__routines
+
+python module zgges__user__routines
@pv Owner
pv added a note

Here you could use f2py's templating:

The syntax <_arg=arg1\,arg2,\0,arg,\2> means the same as <_arg=arg1\,arg2,arg1\,arg2,arg,arg> -- one of the slashes escapes the comma, and apparently the syntax \0, \1, \2, ... can be used (but IMO shouldn't be because it's confusing) to refer to other values in the expansion.

< ftype> :: <_arg>

expands to

real :: arg1,arg2
double precision :: arg1, arg2
complex :: arg
double complex :: arg

I'd maybe also not name the dummy arguments in such a complicated way --- someone reading this code later on will be confused why there are so many underscores, and whether they are needed by f2py for something magical...

@jseabold Collaborator
jseabold added a note

Sorry forgot to look into this. Thanks for the pointer. This is still very murky to me. I still don't see why arg, \2 expands to arg2, arg, arg and why arg, arg should be assumed to be a complex number.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ interface
+ function zselect(e_alpha__i__e,e_beta__i__e)
+ complex*16 :: e_alpha__i__e
+ complex*16 :: e_beta__i__e
+ logical :: zselect
+ end function zselect
+ end interface
+end python module zgges__user__routines
+
+python module dgges__user__routines
+ interface
+ function dselect(e_alphar__i__e,e_alphai__i__e,e_beta__i__e)
+ double precision :: e_alphar__i__e
+ double precision :: e_alphai__i__e
+ double precision :: e_beta__i__e
+ logical :: dselect
+ end function dselect
+ end interface
+end python module dgges__user__routines
+
+python module cgges__user__routines
+ interface
+ function cselect(e_alpha__i__e,e_beta__i__e)
+ complex :: e_alpha__i__e
+ complex :: e_beta__i__e
+ logical :: zselect
+ end function zselect
+ end interface
+end python module cgges__user__routines
+
+python module sgges__user__routines
+ interface
+ function sselect(e_alphar__i__e,e_alphai__i__e,e_beta__i__e)
+ real :: e_alphar__i__e
+ real :: e_alphai__i__e
+ real :: e_beta__i__e
+ logical :: dselect
+ end function dselect
+ end interface
+end python module sgges__user__routines
View
131 scipy/linalg/tests/test_decomp.py
@@ -16,7 +16,8 @@
from scipy.linalg import eig, eigvals, lu, svd, svdvals, cholesky, qr, \
schur, rsf2csf, lu_solve, lu_factor, solve, diagsvd, hessenberg, rq, \
- eig_banded, eigvals_banded, eigh, eigvalsh, qr_multiply, LinAlgError
+ eig_banded, eigvals_banded, eigh, eigvalsh, qr_multiply, LinAlgError, \
+ qz
from scipy.linalg.flapack import dgbtrf, dgbtrs, zgbtrf, zgbtrs, \
dsbev, dsbevd, dsbevx, zhbevd, zhbevx
@@ -1622,7 +1623,133 @@ def test_random_complex(self):
h1 = dot(transp(conj(q)),dot(a,q))
assert_array_almost_equal(h1,h)
-
+class TestQZ(TestCase):
+ def setUp(self):
+ seed(12345)
+
+ def test_qz_single(self):
+ n = 5
+ A = random([n,n]).astype(float32)
+ B = random([n,n]).astype(float32)
+ AA,BB,Q,Z = qz(A,B)
+ assert_array_almost_equal(dot(dot(Q,AA),Z.T), A)
+ assert_array_almost_equal(dot(dot(Q,BB),Z.T), B)
+ assert_array_almost_equal(dot(Q,Q.T), eye(n))
+ assert_array_almost_equal(dot(Z,Z.T), eye(n))
+ assert_(all(diag(BB) >= 0))
+
+ def test_qz_double(self):
+ n = 5
+ A = random([n,n])
+ B = random([n,n])
+ AA,BB,Q,Z = qz(A,B)
+ assert_array_almost_equal(dot(dot(Q,AA),Z.T), A)
+ assert_array_almost_equal(dot(dot(Q,BB),Z.T), B)
+ assert_array_almost_equal(dot(Q,Q.T), eye(n))
+ assert_array_almost_equal(dot(Z,Z.T), eye(n))
+ assert_(all(diag(BB) >= 0))
+
+ def test_qz_complex(self):
+ n = 5
+ A = random([n,n]) + 1j*random([n,n])
+ B = random([n,n]) + 1j*random([n,n])
+ AA,BB,Q,Z = qz(A,B)
+ assert_array_almost_equal(dot(dot(Q,AA),Z.conjugate().T), A)
+ assert_array_almost_equal(dot(dot(Q,BB),Z.conjugate().T), B)
+ assert_array_almost_equal(dot(Q,Q.conjugate().T), eye(n))
+ assert_array_almost_equal(dot(Z,Z.conjugate().T), eye(n))
+ assert_(all(diag(BB) >= 0))
+ assert_(all(diag(BB).imag == 0))
+
+
+ def test_qz_complex64(self):
+ n = 5
+ A = (random([n,n]) + 1j*random([n,n])).astype(complex64)
+ B = (random([n,n]) + 1j*random([n,n])).astype(complex64)
+ AA,BB,Q,Z = qz(A,B)
+ assert_array_almost_equal(dot(dot(Q,AA),Z.conjugate().T), A)
+ assert_array_almost_equal(dot(dot(Q,BB),Z.conjugate().T), B)
+ assert_array_almost_equal(dot(Q,Q.conjugate().T), eye(n))
+ assert_array_almost_equal(dot(Z,Z.conjugate().T), eye(n))
+ assert_(all(diag(BB) >= 0))
+ assert_(all(diag(BB).imag == 0))
+
+ def test_qz_double_complex(self):
+ n = 5
+ A = random([n,n])
+ B = random([n,n])
+ AA,BB,Q,Z = qz(A,B, output='complex')
+ aa = dot(dot(Q,AA),Z.conjugate().T)
+ assert_array_almost_equal(aa.real, A)
+ assert_array_almost_equal(aa.imag, 0)
+ bb = dot(dot(Q,BB),Z.conjugate().T)
+ assert_array_almost_equal(bb.real, B)
+ assert_array_almost_equal(bb.imag, 0)
+ assert_array_almost_equal(dot(Q,Q.conjugate().T), eye(n))
+ assert_array_almost_equal(dot(Z,Z.conjugate().T), eye(n))
+ assert_(all(diag(BB) >= 0))
+
+ def test_qz_double_sort(self):
+ #from http://www.nag.com/lapack-ex/node119.html
+ A = np.array([[3.9, 12.5, -34.5, -0.5],
+ [ 4.3, 21.5, -47.5, 7.5],
+ [ 4.3, 21.5, -43.5, 3.5],
+ [ 4.4, 26.0, -46.0, 6.0 ]])
+
+ B = np.array([[ 1.0, 2.0, -3.0, 1.0],
+ [1.0, 3.0, -5.0, 4.0],
+ [1.0, 3.0, -4.0, 3.0],
+ [1.0, 3.0, -4.0, 4.0]])
+ sort = lambda ar,ai,beta : ai == 0
+
+ AA,BB,Q,Z,sdim = qz(A,B,sort=sort)
+ assert_(sdim == 2)
+ assert_array_almost_equal(dot(dot(Q,AA),Z.T), A)
+ assert_array_almost_equal(dot(dot(Q,BB),Z.T), B)
+
+ # test absolute values bc the sign is ambiguous and might be platform
+ # dependent
+ assert_array_almost_equal(abs(AA), abs(np.array([
+ [3.8009, -69.4505, 50.3135, -43.2884],
+ [0.0000, 9.2033, -0.2001, 5.9881],
+ [0.0000, 0.0000, 1.4279, 4.4453],
+ [0.0000, 0.0000, 0.9019, -1.1962]])), 4)
+ assert_array_almost_equal(abs(BB), abs(np.array([
+ [1.9005, -10.2285, 0.8658, -5.2134],
+ [0.0000, 2.3008, 0.7915, 0.4262],
+ [0.0000, 0.0000, 0.8101, 0.0000],
+ [0.0000, 0.0000, 0.0000, -0.2823]])), 4)
+ assert_array_almost_equal(abs(Q), abs(np.array([
+ [0.4642, 0.7886, 0.2915, -0.2786],
+ [0.5002, -0.5986, 0.5638, -0.2713],
+ [0.5002, 0.0154, -0.0107, 0.8657],
+ [0.5331, -0.1395, -0.7727, -0.3151]])), 4)
+ assert_array_almost_equal(dot(Q,Q.T), eye(4))
+ assert_array_almost_equal(abs(Z), abs(np.array([
+ [0.9961, -0.0014, 0.0887, -0.0026],
+ [0.0057, -0.0404, -0.0938, -0.9948],
+ [0.0626, 0.7194, -0.6908, 0.0363],
+ [0.0626, -0.6934, -0.7114, 0.0956]])), 4)
+ assert_array_almost_equal(dot(Z,Z.T), eye(4))
+
+ def test_qz_complex_sort(self):
+ cA = np.array([
+ [-21.10+22.50*1j, 53.50+-50.50*1j, -34.50+127.50*1j, 7.50+ 0.50*1j],
+ [-0.46+ -7.78*1j, -3.50+-37.50*1j, -15.50+ 58.50*1j,-10.50+ -1.50*1j],
+ [ 4.30+ -5.50*1j, 39.70+-17.10*1j, -68.50+ 12.50*1j, -7.50+ -3.50*1j],
+ [ 5.50+ 4.40*1j, 14.40+ 43.30*1j, -32.50+-46.00*1j,-19.00+-32.50*1j]])
+
+ cB = np.array([
+ [1.00+ -5.00*1j, 1.60+ 1.20*1j,-3.00+ 0.00*1j, 0.00+ -1.00*1j],
+ [0.80+ -0.60*1j, 3.00+ -5.00*1j,-4.00+ 3.00*1j,-2.40+ -3.20*1j],
+ [1.00+ 0.00*1j, 2.40+ 1.80*1j,-4.00+ -5.00*1j, 0.00+ -3.00*1j],
+ [0.00+ 1.00*1j,-1.80+ 2.40*1j, 0.00+ -4.00*1j, 4.00+ -5.00*1j]])
+
+ AAS,BBS,QS,ZS,sdim = qz(cA,cB,sort='lhp')
+
+ eigenvalues = diag(AAS)/diag(BBS)
+ assert_(all(np.real(eigenvalues[:sdim] < 0)))
+ assert_(all(np.real(eigenvalues[sdim:] > 0)))
class TestDatacopied(TestCase):
Something went wrong with that request. Please try again.