diff --git a/CMakeLists.txt b/CMakeLists.txt index a1038751..be2c0718 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,6 +19,7 @@ endif() enable_language(C) enable_language(Fortran) + find_package(PythonLibs REQUIRED) find_package(NumPy REQUIRED) find_package(BLAS REQUIRED) diff --git a/slycot/__init__.py b/slycot/__init__.py index ca429748..11b13d12 100644 --- a/slycot/__init__.py +++ b/slycot/__init__.py @@ -22,8 +22,8 @@ # Identification routines (0/5 wrapped) - # Mathematical routines (3/81 wrapped) - from .math import mc01td, mb05md, mb05nd + # Mathematical routines (6/81 wrapped) + from .math import mc01td, mb03vd, mb03vy, mb03wd, mb05md, mb05nd # Synthesis routines (14/50 wrapped) from .synthesis import sb01bd,sb02md,sb02mt,sb02od,sb03md,sb03od diff --git a/slycot/math.py b/slycot/math.py index 8ad0b020..78808210 100644 --- a/slycot/math.py +++ b/slycot/math.py @@ -70,6 +70,361 @@ def mc01td(dico,dp,p): return out[:-2] +def mb03vd(n, ilo, ihi, A): + """ HQ, Tau = mb03vd(n, ilo, ihi, A) + + To reduce a product of p real general matrices A = A_1*A_2*...*A_p + to upper Hessenberg form, H = H_1*H_2*...*H_p, where H_1 is + upper Hessenberg, and H_2, ..., H_p are upper triangular, by using + orthogonal similarity transformations on A, + + Q_1' * A_1 * Q_2 = H_1, + Q_2' * A_2 * Q_3 = H_2, + ... + Q_p' * A_p * Q_1 = H_p. + + Parameters + ---------- + + n : int + The order of the square matrices A_1, A_2, ..., A_p. + n >= 0. + + ilo, ihi : int + It is assumed that all matrices A_j, j = 2, ..., p, are + already upper triangular in rows and columns [:ilo-1] and + [ihi:n], and A_1 is upper Hessenberg in rows and columns + [:ilo-1] and [ihi:n], with A_1[ilo-1,ilo-2] = 0 (unless + ilo = 1), and A_1[ihi,ihi-1] = 0 (unless ihi = n). + If this is not the case, ilo and ihi should be set to 1 + and n, respectively. + 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. + + A : ndarray + A[:n,:n,:p] must contain the matrices of factors to be reduced; + specifically, A[:,:,j-1] must contain A_j, j = 1, ..., p. + + + Returns + ------- + + HQ : ndarray + 3D array with same shape as A. The upper triangle and the first + subdiagonal of HQ[:n,:n,0] contain the upper Hessenberg + matrix H_1, and the elements below the first subdiagonal, + with the first column of the array Tau represent the + orthogonal matrix Q_1 as a product of elementary + reflectors. See FURTHER COMMENTS. + For j > 1, the upper triangle of HQ[:n,:n,j-1] + contains the upper triangular matrix H_j, and the elements + below the diagonal, with the j-th column of the array TAU + represent the orthogonal matrix Q_j as a product of + elementary reflectors. See FURTHER COMMENTS. + + Tau : ndarray + 2D array with shape (max(1, n-1), p). + The leading n-1 elements in the j-th column contain the + scalar factors of the elementary reflectors used to form + the matrix Q_j, j = 1, ..., p. See FURTHER COMMENTS. + + Raises + ------ + + ValueError : e + e.info contains information about the exact type of exception + + Further Comments + ---------------- + + Each matrix Q_j is represented as a product of (ihi-ilo) + elementary reflectors, + + Q_j = H_j(ilo) H_j(ilo+1) . . . H_j(ihi-1). + + Each H_j(i), i = ilo, ..., ihi-1, has the form + + H_j(i) = I - tau_j * v_j * v_j', + + where tau_j is a real scalar, and v_j is a real vector with + v_j(1:i) = 0, v_j(i+1) = 1 and v_j(ihi+1:n) = 0; v_j(i+2:ihi) + is stored on exit in A_j(i+2:ihi,i), and tau_j in TAU(i,j). + + The contents of A_1 are illustrated by the following example + for n = 7, ilo = 2, and ihi = 6: + + on entry on exit + + ( a a a a a a a ) ( a h h h h h a ) + ( 0 a a a a a a ) ( 0 h h h h h a ) + ( 0 a a a a a a ) ( 0 h h h h h h ) + ( 0 a a a a a a ) ( 0 v2 h h h h h ) + ( 0 a a a a a a ) ( 0 v2 v3 h h h h ) + ( 0 a a a a a a ) ( 0 v2 v3 v4 h h h ) + ( 0 0 0 0 0 0 a ) ( 0 0 0 0 0 0 a ) + + where a denotes an element of the original matrix A_1, h denotes + a modified element of the upper Hessenberg matrix H_1, and vi + denotes an element of the vector defining H_1(i). + + The contents of A_j, j > 1, are illustrated by the following + example for n = 7, ilo = 2, and ihi = 6: + + on entry on exit + + ( a a a a a a a ) ( a h h h h h a ) + ( 0 a a a a a a ) ( 0 h h h h h h ) + ( 0 a a a a a a ) ( 0 v2 h h h h h ) + ( 0 a a a a a a ) ( 0 v2 v3 h h h h ) + ( 0 a a a a a a ) ( 0 v2 v3 v4 h h h ) + ( 0 a a a a a a ) ( 0 v2 v3 v4 v5 h h ) + ( 0 0 0 0 0 0 a ) ( 0 0 0 0 0 0 a ) + + where a denotes an element of the original matrix A_j, h denotes + a modified element of the upper triangular matrix H_j, and vi + denotes an element of the vector defining H_j(i). (The element + (1,2) in A_p is also unchanged for this example.) + + """ + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'p' + hidden, + 'ilo', 'ihi', 'a', + 'lda1' + hidden, 'lda2' + hidden, 'tau', + 'ldtau' + hidden, 'dwork' + hidden, 'info'] + + HQ, Tau, info = _wrapper.mb03vd(n, ilo, ihi, A) + + if info != 0: + e = ValueError( + "Argument '{}' had an illegal value".format(arg_list[-info-1])) + e.info = info + raise e + return (HQ, Tau) + + +def mb03vy(n, ilo, ihi, A, Tau, ldwork=None): + """ Q = mb03vy(n, ilo, ihi, A, Tau, [ldwork]) + + To generate the real orthogonal matrices Q_1, Q_2, ..., Q_p, + which are defined as the product of ihi-ilo elementary reflectors + of order n, as returned by SLICOT Library routine MB03VD: + + Q_j = H_j(ilo) H_j(ilo+1) . . . H_j(ihi-1). + + Parameters + ---------- + + n : int + The order of the matrices Q_1, Q_2, ..., Q_p. N >= 0. + + ilo, ihi : int + The values of the indices ilo and ihi, respectively, used + in the previous call of the SLICOT Library routine MB03VD. + 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. + + A : ndarray + A[:n,:n,j-1] must contain the vectors which define the + elementary reflectors used for reducing A_j, as returned + by SLICOT Library routine MB03VD, j = 1, ..., p. + + Tau : ndarray + The leading N-1 elements in the j-th column must contain + the scalar factors of the elementary reflectors used to + form the matrix Q_j, as returned by SLICOT Library routine + MB03VD. + + ldwork : int, optional + The length of the internal array DWORK. ldwork >= max(1, n). + For optimum performance ldwork should be larger. + + + Returns + ------- + + Q : ndarray + 3D array with same shape as A. Q[:n,:n,j-1] contains the + N-by-N orthogonal matrix Q_j, j = 1, ..., p. + + Raises + ------ + + ValueError : + e.info contains the number of the argument that was invalid + + """ + + hidden = ' (hidden by the wrapper)' + arg_list = ['n', 'p' + hidden, + 'ilo', 'ihi', 'a', + 'lda1' + hidden, 'lda2' + hidden, 'tau', + 'ldtau' + hidden, 'dwork' + hidden, 'ldwork', 'info' + hidden] + + if not ldwork: + ldwork = max(1, 2 * n) + + Q, info = _wrapper.mb03vy(n, ilo, ihi, A, Tau, ldwork) + + if info != 0: + e = ValueError( + "Argument '{}' had an illegal value".format(arg_list[-info-1])) + e.info = info + raise e + + return Q + + +def mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork=None): + """ T, Z, Wr = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, H, Q, [ldwork]) + + To compute the Schur decomposition and the eigenvalues of a + product of matrices, H = H_1*H_2*...*H_p, with H_1 an upper + Hessenberg matrix and H_2, ..., H_p upper triangular matrices, + without evaluating the product. Specifically, the matrices Z_i + are computed, such that + + Z_1' * H_1 * Z_2 = T_1, + Z_2' * H_2 * Z_3 = T_2, + ... + Z_p' * H_p * Z_1 = T_p, + + where T_1 is in real Schur form, and T_2, ..., T_p are upper + triangular. + + The routine works primarily with the Hessenberg and triangular + submatrices in rows and columns ilo to ihi, but optionally applies + the transformations to all the rows and columns of the matrices + H_i, i = 1,...,p. The transformations can be optionally + accumulated. + + Parameters + ---------- + + job : {'E', 'S'} + Indicates whether the user wishes to compute the full + Schur form or the eigenvalues only, as follows: + = 'E': Compute the eigenvalues only; + = 'S': Compute the factors T_1, ..., T_p of the full + Schur form, T = T_1*T_2*...*T_p. + + compz : {'N', 'I', 'V'} + Indicates whether or not the user wishes to accumulate + the matrices Z_1, ..., Z_p, as follows: + = 'N': The matrices Z_1, ..., Z_p are not required; + = 'I': Z_i is initialized to the unit matrix and the + orthogonal transformation matrix Z_i is returned, + i = 1, ..., p; + = 'V': Z_i must contain an orthogonal matrix Q_i on + entry, and the product Q_i*Z_i is returned, + i = 1, ..., p. + + n : int + The order of the matrix H. n >= 0 + + ilo, ihi : int + It is assumed that all matrices H_j, j = 2, ..., p, are + already upper triangular in rows and columns [:ilo-1] and + [ihi:n], and H_1 is upper quasi-triangular in rows and + columns [:ilo-1] and [ihi:n], with H_1[ilo-1,ilo-2] = 0 + (unless ilo = 1), and H_1[ihi,ihi-1] = 0 (unless ihi = n). + The routine works primarily with the Hessenberg submatrix + in rows and columns ilo to ihi, but applies the + transformations to all the rows and columns of the + matrices H_i, i = 1,...,p, if JOB = 'S'. + 1 <= ilo <= max(1,n); min(ilo,n) <= ihi <= n. + + iloz, ihiz : int + Specify the rows of Z to which the transformations must be + applied if compz = 'I' or compz = 'V'. + 1 <= iloz <= ilo; ihi <= ihiz <= n. + + H : ndarray + H[:n,:n,0] must contain the upper Hessenberg matrix H_1 and + H[:n,:n,j-1] for j > 1 must contain the upper triangular matrix + H_j, j = 2, ..., p. + + Q : ndarray + If compz = 'V', Q[:n,:n,:p] must contain the current matrix Q of + transformations accumulated by SLICOT Library routine + MB03VY. + If compz = 'I', Q is ignored + + ldwork : int, optinal + The length of the cache array. The default value is + ihi-ilo+p-1 + + + + Returns + ------- + + T : ndarray + 3D array with the same shape as H. + If JOB = 'S', T[:n,:n,0] is upper quasi-triangular in rows + and columns [ilo-1:ihi], with any 2-by-2 diagonal blocks + corresponding to a pair of complex conjugated eigenvalues, and + T[:n,:n,j-1] for j > 1 contains the resulting upper + triangular matrix T_j. + If job = 'E', T is None + + Z : ndarray + 3D array with the same shape as Q. + If compz = 'V', or compz = 'I', the leading + N-by-N-by-P part of this array contains the transformation + matrices which produced the Schur form; the + transformations are applied only to the submatrices + Z[iloz-1:ihiz,ilo-1:ihi,j-1], j = 1, ..., p. + If compz = 'N', Z is None + + W : ndarray (dtype=complex) + 1D array with shape (n). + The computed eigenvalues ilo to ihi. If two eigenvalues + are computed as a complex conjugate pair, they are stored + in consecutive elements of W say the i-th and + (i+1)th, with imag(W][i]) > 0 and imag(W[i+1]) < 0. + If JOB = 'S', the eigenvalues are stored in the same order + as on the diagonal of the Schur form returned in H. + + Raises + ------ + + ValueError : e + e.info contains information about the exact type of exception + + """ + hidden = ' (hidden by the wrapper)' + arg_list = ['job', 'compz', 'n', 'p' + hidden, + 'ilo', 'ihi', 'iloz', 'ihiz', + 'h', 'ldh1' + hidden, 'ldh2' + hidden, + 'z', 'ldz1' + hidden, 'ldz2' + hidden, + 'wr', 'wi', + 'dwork' + hidden, 'ldwork', 'info' + hidden] + + if not ldwork: + p = H.shape[2] + ldwork = max(1, ihi - ilo + p - 1) + + T, Z, Wr, Wi, info = _wrapper.mb03wd( + job, compz, n, ilo, ihi, iloz, ihiz, H, Q, ldwork) + + if info < 0: + e = ValueError( + "Argument '{}' had an illegal value".format(arg_list[-info-1])) + e.info = info + raise e + elif info > 0: + warnings.warn(("failed to compute all the eigenvalues {ilo} to {ihi} " + "in a total of 30*({ihi}-{ilo}+1) iterations " + "the elements {i}:{ihi} of Wr contain those " + "eigenvalues which have been successfully computed." + ).format(i=info, ilo=ilo, ihi=ihi)) + if job == 'E': + T = None + if compz == 'N': + Z = None + + W = Wr + Wi*1J + return (T, Z, W) + + def mb05md(a, delta, balanc='N'): """Ar, Vr, Yr, VAL = mb05md(a, delta, balanc='N') @@ -166,7 +521,7 @@ def mb05md(a, delta, balanc='N'): from slycot import mb05nd import numpy as np a = np.mat('[-2. 0; 0.1 -3.]') -mb05nd(a.shape[0], a, 0.1) +mb05nd(a.shape[0], a, 0.1) """ def mb05nd(a, delta, tol=1e-7): diff --git a/slycot/src/math.pyf b/slycot/src/math.pyf index 101e20e2..ecf94f2c 100644 --- a/slycot/src/math.pyf +++ b/slycot/src/math.pyf @@ -12,6 +12,57 @@ subroutine mc01td(dico,dp,p,stable,nz,dwork,iwarn,info) ! in :new:MC01TD.f integer intent(out) :: info end subroutine mc01td +subroutine mb03vd(n,p,ilo,ihi,a,lda1,lda2,tau,ldtau,dwork,info) ! in MB03VD.f + integer intent(in),check(n>=0) :: n + integer intent(hide),depend(a),check(p>=1) :: p=shape(a,2) + integer intent(in),depend(n),check(1<=ilo && ilo<=max(1,n)) :: ilo + integer intent(in),depend(n,ilo),check(min(ilo,n)<=ihi && ihi<=n) :: ihi + double precision intent(in,out,copy),dimension(lda1,lda2,p) :: a + integer intent(hide),depend(a,n),check(lda1>=max(1,n)) :: lda1=shape(a,0) + integer intent(hide),depend(a,n),check(lda2>=max(1,n)) :: lda2=shape(a,1) + double precision intent(out),depend(n,p),dimension(max(1,n-1),p) :: tau + integer intent(hide),depend(tau) :: ldtau=shape(tau,0) + double precision intent(hide,cache),depend(n),dimension(n) :: dwork + integer intent(out) :: info +end subroutine mb03vd + +subroutine mb03vy(n,p,ilo,ihi,a,lda1,lda2,tau,ldtau,dwork,ldwork,info) ! in MB03VY.f + integer intent(in),check(n>=0) :: n + integer intent(hide),depend(a),check(p>=1) :: p=shape(a,2) + integer intent(in),depend(n),check(1<=ilo && ilo<=max(1,n)) :: ilo + integer intent(in),depend(n,ilo),check(min(ilo,n)<=ihi && ihi<=n) :: ihi + double precision intent(in,out,copy),dimension(lda1,lda2,p) :: a + integer intent(hide),depend(a,n),check(lda1>=max(1,n)) :: lda1=shape(a,0) + integer intent(hide),depend(a,n),check(lda2>=max(1,n)) :: lda2=shape(a,1) + double precision intent(in),depend(n,p),dimension(ldtau,p) :: tau + integer intent(hide),depend(tau),check(ldtau>=max(1,n-1)) :: ldtau=shape(tau,0) + double precision intent(hide,cache),depend(ldwork),dimension(ldwork) :: dwork + integer intent(in),optional,check(ldwork>=max(1,n)) :: ldwork=max(1,n) + integer intent(out) :: info +end subroutine mb03vy + +subroutine mb03wd(job,compz,n,p,ilo,ihi,iloz,ihiz,h,ldh1,ldh2,z,ldz1,ldz2,wr,wi,dwork,ldwork,info) ! in MB03WD.f + character intent(in) :: job + character intent(in) :: compz + integer intent(in),check(n>=0) :: n + integer intent(hide),depend(h),check(p>=1) :: p=shape(h,2) + integer intent(in),depend(n),check(1<=ilo && ilo<=max(1,n)) :: ilo + integer intent(in),depend(n,ilo),check(min(ilo,n)<=ihi && ihi<=n) :: ihi + integer intent(in),depend(ilo),check(1<=iloz & iloz<=ilo) :: iloz + integer intent(in),depend(n,ihi),check(ihi<=ihiz && ihiz<=n) :: ihiz + double precision intent(in,out,copy),dimension(ldh1,ldh2,p) :: h + integer intent(hide),depend(h,n),check(ldh1>=max(1,n)) :: ldh1=shape(h,0) + integer intent(hide),depend(h,n),check(ldh2>=max(1,n)) :: ldh2=shape(h,1) + double precision intent(in,out,copy),depend(p),dimension(ldz1,ldz2,p) :: z + integer intent(hide),depend(z) :: ldz1=shape(z,0) + integer intent(hide),depend(z) :: ldz2=shape(z,1) + double precision intent(out), dimension(n), depend(n) :: wr + double precision intent(out), dimension(n), depend(n) :: wi + double precision intent(hide,cache), dimension(ldwork) :: dwork + integer optional,check(ldwork>=ihi-ilo+p-1), depend(ihi,ilo,p) :: ldwork=max(1,ihi-ilo+p-1) + integer intent(out) :: info +end subroutine mb03wd + subroutine mb05md(balanc,n,delta,a,lda,v,ldv,y,ldy,valr,vali,iwork,dwork,ldwork,info) ! in MB05MD.f character intent(in):: balanc integer intent(in),check(n>=0) :: n @@ -45,3 +96,6 @@ subroutine mb05nd(n,delta,a,lda,ex,ldex,exint,ldexin,tol,iwork,dwork,ldwork,info integer intent(hide),depend(dwork):: ldwork=shape(dwork,0) integer intent(out) :: info end subroutine mb05nd + +! This file was auto-generated with f2py (version:2). +! See http://cens.ioc.ee/projects/f2py2e/ diff --git a/slycot/tests/test_mb.py b/slycot/tests/test_mb.py index 2f8e3750..b55a0d91 100755 --- a/slycot/tests/test_mb.py +++ b/slycot/tests/test_mb.py @@ -6,13 +6,152 @@ import unittest import numpy as np -from slycot import mb05md, mb05nd +from slycot import mb03vd, mb03vy, mb03wd, mb05md, mb05nd from numpy.testing import assert_allclose class test_mb(unittest.TestCase): + def test_mb03vd_mb03vy_ex(self): + """Test MB03VD and MB03VY + with the example given in the MB03VD SLICOT documentation""" + + n = 4 + p = 2 + ilo = 1 + ihi = 4 + A = np.zeros((n, n, p)) + A[:, :, 0] = [[1.5, -.7, 3.5, -.7], + [1. , 0. , 2. , 3. ], + [1.5, -.7, 2.5, -.3], + [1. , 0. , 2. , 1. ]] + A[:, :, 1] = [[1.5, -.7, 3.5, -.7], + [1. , 0. , 2. , 3. ], + [1.5, -.7, 2.5, -.3], + [1. , 0. , 2. , 1. ]] + + H_ref = np.zeros((n, n, p)) + H_ref[:, :, 0] = [[-2.3926, 2.7042, -0.9598, -1.2335], + [ 4.1417, -1.7046, 1.3001, -1.3120], + [ 0.0000, -1.6247, -0.2534, 1.6453], + [ 0.0000, 0.0000, -0.0169, -0.4451]] + + H_ref[:, :, 1] = [[-2.5495, 2.3402, 4.7021, 0.2329], + [ 0.0000, 1.9725, -0.2483, -2.3493], + [ 0.0000, 0.0000, -0.6290, -0.5975], + [ 0.0000, 0.0000, 0.0000, -0.4426]] + + Q_ref = np.zeros((n, n, p)) + Q_ref[:, :, 0] = [[ 1.0000, 0.0000, 0.0000, 0.0000], + [ 0.0000, -0.7103, 0.5504, -0.4388], + [ 0.0000, -0.4735, -0.8349, -0.2807], + [ 0.0000, -0.5209, 0.0084, 0.8536]] + + Q_ref[:, :, 1] = [[-0.5883, 0.2947, 0.7528, -0.0145], + [-0.3922, -0.8070, 0.0009, -0.4415], + [-0.5883, 0.4292, -0.6329, -0.2630], + [-0.3922, -0.2788, -0.1809, 0.8577]] + + HQ, Tau = mb03vd(n, ilo, ihi, A) + + H = np.zeros_like(HQ) + Q = np.zeros_like(HQ) + + for k in range(p): + Q[:, :, k] = np.tril(HQ[:, :, k]) + if k == 0: + H[:, :, k] = np.triu(HQ[:n, :n, k], -1) + elif k > 0: + H[:, :, k] = np.triu(HQ[:n, :n, k]) + assert_allclose(H[:, :, k], H_ref[:, :, k], atol=1e-4) + + Qr = mb03vy(n, ilo, ihi, Q, Tau) + + for k in range(p): + assert_allclose(Qr[:, :, k], Q_ref[:, :, k], atol=1e-4) + + # Computer Error: too machine dependent to test to reference value + # SSQ_ref = 2.93760e-15 + # SSQ = 0. + # for k in range(p): + # kp1 = k+1 + # if kp1 > p-1: + # kp1 = 0 + # P = Qr[:, :, k].T.dot(A[: ,: ,k]).dot(Qr[: ,: ,kp1]) - H[: ,: ,k] + # SSQ = np.sqrt(SSQ**2 + np.linalg.norm(P,'fro')**2) + + def test_mb03wd_ex(self): + """Test MB03WD with the example given in the SLICOT documentation""" + + n = 4 + p = 2 + ilo = 1 + ihi = 4 + iloz = 1 + ihiz = 4 + job = 'S' + compz = 'V' + A = np.zeros((n, n, p)) + A[:, :, 0] = [[1.5, -.7, 3.5, -.7], + [1. , 0. , 2. , 3. ], + [1.5, -.7, 2.5, -.3], + [1. , 0. , 2. , 1. ]] + A[:, :, 1] = [[1.5, -.7, 3.5, -.7], + [1. , 0. , 2. , 3. ], + [1.5, -.7, 2.5, -.3], + [1. , 0. , 2. , 1. ]] + + W_ref = np.array([6.449861+7.817717J, + 6.449861-7.817717J, + 0.091315+0.000000J, + 0.208964+0.000000J]) + + T_ref = np.zeros((n, n, p)) + T_ref[:, :, 0] = [[ 2.2112, 4.3718, -2.3362, 0.8907], + [ -0.9179, 2.7688, -0.6570, -2.2426], + [ 0.0000, 0.0000, 0.3022, 0.1932], + [ 0.0000, 0.0000, 0.0000, -0.4571]] + + T_ref[:, :, 1] = [[ 2.9169, 3.4539, 2.2016, 1.2367], + [ 0.0000, 3.4745, 1.0209, -2.0720], + [ 0.0000, 0.0000, 0.3022, -0.1932], + [ 0.0000, 0.0000, 0.0000, -0.4571]] + + Z_ref = np.zeros((n, n, p)) + Z_ref[:, :, 0] = [[ 0.3493, 0.6751, -0.6490, 0.0327], + [ 0.7483, -0.4863, -0.1249, -0.4336], + [ 0.2939, 0.5504, 0.7148, -0.3158], + [ 0.4813, -0.0700, 0.2286, 0.8433]] + + + Z_ref[:, :, 1] = [[ 0.2372, 0.7221, 0.6490, 0.0327], + [ 0.8163, -0.3608, 0.1249, -0.4336], + [ 0.2025, 0.5902, -0.7148, -0.3158], + [ 0.4863, 0.0076, -0.2286, 0.8433]] + + HQ, Tau = mb03vd(n, ilo, ihi, A) + Q = mb03vy(n, ilo, ihi, HQ, Tau) + T, Z, W = mb03wd(job, compz, n, ilo, ihi, iloz, ihiz, HQ, Q) + + # TODO (?) + # isolate eigenvalues with math.mb03wx + + assert_allclose(W, W_ref, atol=1e-5) + assert_allclose(T, T_ref, atol=1e-4) + assert_allclose(Z, Z_ref, atol=1e-4) + + # Computer Error: too machine dependent to test to reference value + # SSQ_ref = 7.18432D-15 + # SSQ = 0. + # for k in range(p): + # kp1 = k+1 + # if kp1 > p-1: + # kp1 = 0 + # P = Zrr[:, :, k].T.dot(A[: ,: ,k]).dot(Zrr[: ,: ,kp1]) - Hrr[: ,: ,k] + # SSQ = np.sqrt(SSQ**2 + np.linalg.norm(P,'fro')**2) + + def test_mb05md(self): """ test_mb05md: verify Matrix exponential with slicot doc example data from http://slicot.org/objects/software/shared/doc/MB05MD.html