adding the mb03rd wrapper

repagh · repagh · commit 009a2b804ee7 · 2019-08-02T22:34:28.000+02:00
diff --git a/slycot/src/transform.pyf b/slycot/src/transform.pyf
@@ -528,3 +528,40 @@ subroutine tg01fd_uu(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ld
 integer required intent(in) :: ldwork
     integer intent(out) :: info
 end subroutine tg01fd_uu
+subroutine mb03rd_n(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
+  fortranname mb03rd
+  character intent(hide) :: jobx = 'N'
+  character intent(in),required :: sort
+  integer intent(in),required,check(n>0) :: n
+  double precision intent(in),required,check(pmax>=1.0) :: pmax
+  double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
+  integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
+  double precision intent(cache,hide) :: x
+  integer intent(in,hide) :: ldx=1
+  integer intent(out) :: nblcks
+  integer intent(out),dimension(n),depend(n) :: blsize
+  double precision intent(out),dimension(n),depend(n) :: wr
+  double precision intent(out),dimension(n),depend(n) :: wi
+  double precision intent(in) :: tol
+  double precision intent(cache,hide),dimension(n),depend(n) :: dwork
+  integer intent(out) :: info
+end subroutine mb03rd_n
+subroutine mb03rd_u(jobx,sort,n,pmax,a,lda,x,ldx,nblcks,blsize,wr,wi,tol,dwork,info) ! in MB03RD.f
+  fortranname mb03rd
+  character intent(hide) :: jobx = 'U'
+  character intent(in),required :: sort
+  integer intent(in),required,check(n>0) :: n
+  double precision intent(in),required,check(pmax>=1.0) :: pmax
+  double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
+  integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
+  double precision intent(in,out,copy),dimension(n,n),depend(n) :: x
+  integer intent(hide),depend(x) :: ldx=MAX(shape(x,0),1)
+  integer intent(out) :: nblcks
+  integer intent(out),dimension(n),depend(n) :: blsize
+  double precision intent(out),dimension(n),depend(n) :: wr
+  double precision intent(out),dimension(n),depend(n) :: wi
+  double precision intent(in) :: tol
+  double precision intent(cache,hide),dimension(n),depend(n) :: dwork
+  integer intent(out) :: info
+end subroutine mb03rd_u
+
diff --git a/slycot/tests/test_mb03rd.py b/slycot/tests/test_mb03rd.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+#
+# test_mb03rd.py - test suite for Shur form reduction
+# RvP, 31 Jul 2019
+import unittest
+from slycot import transform
+import numpy as np
+from numpy.testing import assert_raises, assert_almost_equal, assert_equal
+from scipy.linalg import schur
+
+test1_A = np.array([
+   [ 1.,  -1.,   1.,   2.,   3.,   1.,   2.,   3.],
+   [ 1.,   1.,   3.,   4.,   2.,   3.,   4.,   2.],
+   [ 0.,   0.,   1.,  -1.,   1.,   5.,   4.,   1.],
+   [ 0.,   0.,   0.,   1.,  -1.,   3.,   1.,   2.],
+   [ 0.,   0.,   0.,   1.,   1.,   2.,   3.,  -1.],
+   [ 0.,   0.,   0.,   0.,   0.,   1.,   5.,   1.],
+   [ 0.,   0.,   0.,   0.,   0.,   0.,   0.99999999, -0.99999999 ],
+   [ 0.,   0.,   0.,   0.,   0.,   0.,   0.99999999,  0.99999999 ]
+   ])
+test1_n = test1_A.shape[0]
+
+test1_Ar = np.array([
+    [ 1.0000, -1.0000, -1.2247, -0.7071, -3.4186,  1.4577,  0.0000,  0.0000 ],
+    [ 1.0000,  1.0000,  0.0000,  1.4142, -5.1390,  3.1637,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  1.0000, -1.7321, -0.0016,  2.0701,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  0.5774,  1.0000,  0.7516,  1.1379,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  1.0000, -5.8606,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  0.1706,  1.0000,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  1.0000, -0.8850 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  1.0000 ],
+    ])
+
+test1_Xr = np.array([
+    [ 1.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.9045,  0.1957 ],
+    [ 0.0000,  1.0000,  0.0000,  0.0000,  0.0000,  0.0000, -0.3015,  0.9755 ],
+    [ 0.0000,  0.0000,  0.8165,  0.0000, -0.5768, -0.0156, -0.3015,  0.0148 ],
+    [ 0.0000,  0.0000, -0.4082,  0.7071, -0.5768, -0.0156,  0.0000, -0.0534 ],
+    [ 0.0000,  0.0000, -0.4082, -0.7071, -0.5768, -0.0156,  0.0000,  0.0801 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000, -0.0276,  0.9805,  0.0000,  0.0267 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0332, -0.0066,  0.0000,  0.0000 ],
+    [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0011,  0.1948,  0.0000,  0.0000 ]
+    ])
+
+test1_pmax = 1e3
+test1_tol = 0.01
+class test_mb03rd(unittest.TestCase):
+    def test1(self):
+        # create schur form with scipy
+        A, X = schur(test1_A)
+        Ah, Xh = np.copy(A), np.copy(X)
+        # on this basis, get the transform
+        Ar, Xr, blks, eig = transform.mb03rd(
+            test1_n, A, X, 'U', 'S', test1_pmax, test1_tol)
+        # ensure X and A are unchanged
+        assert_equal(A, Ah)
+        assert_equal(X, Xh)
+        # compare to test case results
+        assert_almost_equal(Ar, test1_Ar, decimal=4)
+        assert_almost_equal(Xr, test1_Xr, decimal=4)
+        
+if __name__ == "__main__":
+    unittest.main()
diff --git a/slycot/transform.py b/slycot/transform.py
@@ -1353,4 +1353,110 @@ def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ld
     
     return A,E,B,C,ranke,rnka22,Q,Z
 
+def mb03rd(n,A,X=None,jobx='U',sort='N',pmax=1.0,tol=0.0):
+    """ A,X,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol]) -- if jobx='U'
+        A,blcks,EIG = mb03rd(n,A,[X,job,sort,pmax,tol])   -- if jobx='N'
+
+    To reduce a matrix A in real Schur form to a block-diagonal form
+    using well-conditioned non-orthogonal similarity transformations.
+    The condition numbers of the transformations used for reduction
+    are roughly bounded by pmax*pmax, where pmax is a given value.
+    The transformations are optionally postmultiplied in a given
+    matrix X. The real Schur form is optionally ordered, so that
+    clustered eigenvalues are grouped in the same block.
+
+    Required arguments:
+        n : input int
+            The order of the matrices A and X.  n >= 0.
+        A : input rank-2 array('d') with bounds (n,n)
+            the matrix A to be block-diagonalized, in real Schur form.
+    Optional arguments:
+        X : input rank-2 array('d') with bounds (n,n)
+            a given matrix X, for accumulation of transformations (only if
+            jobx='U'
+        jobx : input char*1
+            Specifies whether or not the transformations are
+            accumulated, as follows:
+            = 'N':  The transformations are not accumulated;
+            = 'U':  The transformations are accumulated in X (the
+                    given matrix X is updated).
+        sort : input char*1
+            Specifies whether or not the diagonal blocks of the real
+            Schur form are reordered, as follows:
+            = 'N':  The diagonal blocks are not reordered;
+            = 'S':  The diagonal blocks are reordered before each
+                    step of reduction, so that clustered eigenvalues
+                    appear in the same block;
+            = 'C':  The diagonal blocks are not reordered, but the
+                    "closest-neighbour" strategy is used instead of
+                    the standard "closest to the mean" strategy
+                    (see METHOD);
+            = 'B':  The diagonal blocks are reordered before each
+                    step of reduction, and the "closest-neighbour"
+                    strategy is used (see METHOD).
+        pmax : input float
+            An upper bound for the infinity norm of elementary
+            submatrices of the individual transformations used for
+            reduction (see METHOD).  PMAX >= 1.0D0.
+        tol : input float
+            The tolerance to be used in the ordering of the diagonal
+            blocks of the real Schur form matrix.
+            If the user sets TOL > 0, then the given value of TOL is
+            used as an absolute tolerance: a block i and a temporarily
+            fixed block 1 (the first block of the current trailing
+            submatrix to be reduced) are considered to belong to the
+            same cluster if their eigenvalues satisfy
+
+              | lambda_1 - lambda_i | <= TOL.
+
+            If the user sets TOL < 0, then the given value of TOL is
+            used as a relative tolerance: a block i and a temporarily
+            fixed block 1 are considered to belong to the same cluster
+            if their eigenvalues satisfy, for j = 1, ..., N,
+
+              | lambda_1 - lambda_i | <= | TOL | * max | lambda_j |.
+
+            If the user sets TOL = 0, then an implicitly computed,
+            default tolerance, defined by TOL = SQRT( SQRT( EPS ) )
+            is used instead, as a relative tolerance, where EPS is
+            the machine precision (see LAPACK Library routine DLAMCH).
+            If SORT = 'N' or 'C', this parameter is not referenced.
+    Return objects:
+        Ar : output rank-2 array('d') with bounds (n,n)
+            Contains the computed block-diagonal matrix, in real Schur
+            canonical form. The non-diagonal blocks are set to zero.
+        Xr : output rank-2 array('d') with bounds (n,n)
+            Contains the product of the given matrix X and the
+            transformation matrix that reduced A to block-diagonal
+            form. The transformation matrix is itself a product of
+            non-orthogonal similarity transformations having elements
+            with magnitude less than or equal to PMAX.
+            If JOBX = 'N', this array is not referenced, and not returned
+        blksize : output rank-1 array('i') with bounds (n)
+            The orders of the resulting diagonal blocks of the matrix Ar.
+        W : output rank-1 array('c') size (n)
+            This arrays contain the eigenvalues of the matrix A.
+"""
+    hidden = ' (hidden by the wrapper)'
+    arg_list = ('jobx', 'sort', 'n', 'pmax', 'A', 'LDA'+hidden,
+                'X', 'LDX'+hidden, 'nblks'+hidden, 'blsize'+hidden,
+                'WR'+hidden, 'WI'+hidden, 'tol',
+                'DWORK'+hidden, 'INFO'+hidden)
+    if jobx == 'N':
+        out = _wrapper.mb03rd_n(sort, n, pmax, A, tol)
+    else:
+        if X is None:
+            X = _np.eye(n)
+        out = _wrapper.mb03rd_u(sort, n, pmax, A, X, tol)
+
+    if out[-1] < 0:
+        error_text = "The following argument had an illegal value: "+arg_list[-out[-1]-1]
+        e = ValueError(error_text)
+        e.info = out[-1]
+        raise e
+    if jobx == 'N':
+        return out[0], out[2][:out[1]], out[-3] + out[-2]*1j
+    else:
+        return out[0], out[1], out[3][:out[2]], out[-3] + out[-2]*1j
+    
 # to be replaced by python wrappers
