Skip to content

Commit 7f5dc1c

Browse files
murrayrmmurrayrm
authored
Merge pull request #36 from marcuslil/master
add support for tg01ad
2 parents cd724f0 + 3de6783 commit 7f5dc1c

File tree

3 files changed

+217
-0
lines changed

3 files changed

+217
-0
lines changed

‎slycot/src/transform.pyf‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -421,6 +421,26 @@ subroutine tb01pd(job,equil,n,m,p,a,lda,b,ldb,c,ldc,nr,tol,iwork,dwork,ldwork,in
421421
integer optional,depend(n,m,p) :: ldwork=max(1,n+max(n,max(3*m,3*p)))
422422
integer intent(out) :: info
423423
end subroutine tb01pd
424+
subroutine tg01ad(job,l,n,m,p,thresh,a,lda,e,lde,b,ldb,c,ldc,lscale,rscale,dwork,info) ! in TG01AD.f
425+
character :: job
426+
integer intent(in),required,check(l>=0) :: l
427+
integer intent(in),required,check(n>=0) :: n
428+
integer intent(in),required,check(m>=0) :: m
429+
integer intent(in),required,check(p>=0) :: p
430+
double precision intent(in),required,check(thresh>=0) :: thresh
431+
double precision intent(in,out,copy),dimension(l,n),depend(l,n) :: a
432+
integer intent(hide),depend(a) :: lda=MAX(shape(a,0),1)
433+
double precision intent(in,out,copy),dimension(l,n),depend(l,m) :: e
434+
integer intent(hide),depend(e) :: lde=MAX(shape(e,0),1)
435+
double precision intent(in,out,copy),dimension(l,m),depend(l,m) :: b
436+
integer intent(hide),depend(b) :: ldb=MAX(shape(b,0),1)
437+
double precision intent(in,out,copy),dimension(p,n),depend(p,n) :: c
438+
integer intent(hide),depend(c) :: ldc=MAX(shape(c,0),1)
439+
double precision intent(out),dimension(l) :: lscale
440+
double precision intent(out),dimension(n) :: rscale
441+
double precision intent(cache,hide),dimension(3*(l+n)) :: dwork
442+
integer intent(out) :: info
443+
end subroutine tg01ad
424444
subroutine tg01fd_nn(compq,compz,joba,l,n,m,p,a,lda,e,lde,b,ldb,c,ldc,q,ldq,z,ldz,ranke,rnka22,tol,iwork,dwork,ldwork,info) ! in TG01FD.f
425445
fortranname tg01fd
426446
character intent(hide) :: compq = 'N'

‎slycot/tests/test_tg01ad.py‎

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
# ===================================================
2+
# tg01ad tests
3+
4+
import unittest
5+
from slycot import transform
6+
import numpy as np
7+
8+
from numpy.testing import assert_raises, assert_almost_equal, assert_equal
9+
10+
# test1 input parameters
11+
12+
test1_l = 4
13+
test1_n = 4
14+
test1_m = 2
15+
test1_p = 2
16+
test1_job = 'A'
17+
test1_thresh = 0.0
18+
19+
test1_A = \
20+
np.array([[-1.0, 0.0, 0.0, 3e-3 ],
21+
[ 0.0, 0.0, 0.1, 2e-2 ],
22+
[ 1e2, 10.0, 0.0, 0.4 ],
23+
[ 0.0, 0.0, 0.0, 0.0 ]])
24+
25+
test1_E = \
26+
np.array([[ 1.0, 0.2, 0.0, 0.0 ],
27+
[ 0.0, 1.0, 0.0, 1e-2 ],
28+
[ 3e2, 90.0, 6.0, 0.3 ],
29+
[ 0.0, 0.0, 20.0, 0.0 ]])
30+
31+
test1_B = \
32+
np.array([[ 10.0, 0.0 ],
33+
[ 0.0, 0.0 ],
34+
[ 0.0, 1e3 ],
35+
[ 1e4, 1e4 ]])
36+
37+
test1_C = \
38+
np.array([[-0.1, 0.0, 1e-3, 0.0 ],
39+
[ 0.0, 1e-2, -1e-3, 1e-4 ]])
40+
41+
test1_A_desired = \
42+
np.array([[-1.0, 0.0, 0.0, 0.3 ],
43+
[ 0.0, 0.0, 1.0, 2.0 ],
44+
[ 1.0, 0.1, 0.0, 0.4 ],
45+
[ 0.0, 0.0, 0.0, 0.0 ]])
46+
47+
test1_E_desired = \
48+
np.array([[ 1.0, 0.2, 0.0, 0.0 ],
49+
[ 0.0, 1.0, 0.0, 1.0 ],
50+
[ 3.0, 0.9, 0.6, 0.3 ],
51+
[ 0.0, 0.0, 0.2, 0.0 ]])
52+
53+
test1_B_desired = \
54+
np.array([[ 1e2, 0.0 ],
55+
[ 0.0, 0.0 ],
56+
[ 0.0, 1e2 ],
57+
[ 1e2, 1e2 ]])
58+
59+
test1_C_desired = \
60+
np.array([[-1e-2, 0.0, 1e-3, 0.0 ],
61+
[ 0.0, 1e-3, -1e-3, 1e-3 ]])
62+
63+
test1_lscale_desired = \
64+
np.array([ 10.0, 10.0, 0.1, 1e-2 ])
65+
66+
test1_rscale_desired = \
67+
np.array([ 0.1, 0.1, 1.0, 10.0 ])
68+
69+
class test_tg01ad(unittest.TestCase):
70+
""" test1: Verify tg01ad with input parameters according to example in documentation """
71+
72+
def test1_tg01ad(self):
73+
74+
A,E,B,C,lscale,rscale = transform.tg01ad(l=test1_l,n=test1_n,m=test1_m,p=test1_p,A=test1_A,E=test1_E,B=test1_B,C=test1_C,job=test1_job, thresh=test1_thresh)
75+
76+
assert_almost_equal(A, test1_A_desired)
77+
assert_almost_equal(E, test1_E_desired)
78+
assert_almost_equal(B, test1_B_desired)
79+
assert_almost_equal(C, test1_C_desired)
80+
assert_almost_equal(lscale, test1_lscale_desired)
81+
assert_almost_equal(rscale, test1_rscale_desired)
82+
83+
84+
def suite():
85+
return unittest.TestLoader().loadTestsFromTestCase(TestConvert)
86+
87+
88+
if __name__ == "__main__":
89+
unittest.main()

‎slycot/transform.py‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1050,6 +1050,114 @@ def tb01pd(n, m, p, A, B, C, job='M', equil='S', tol=1e-8, ldwork=None):
10501050
raise e
10511051
return out[:-1]
10521052

1053+
def tg01ad(l,n,m,p,A,E,B,C,thresh=0.0,job='A'):
1054+
""" A,E,B,C,lscale,rscale = tg01ad(l,n,m,p,A,E,B,C,[thresh,job])
1055+
1056+
To balance the matrices of the system pencil
1057+
1058+
S = ( A B ) - lambda ( E 0 ) := Q - lambda Z,
1059+
( C 0 ) ( 0 0 )
1060+
1061+
corresponding to the descriptor triple (A-lambda E,B,C),
1062+
by balancing. This involves diagonal similarity transformations
1063+
(Dl*A*Dr - lambda Dl*E*Dr, Dl*B, C*Dr) applied to the system
1064+
(A-lambda E,B,C) to make the rows and columns of system pencil
1065+
matrices
1066+
1067+
diag(Dl,I) * S * diag(Dr,I)
1068+
1069+
as close in norm as possible. Balancing may reduce the 1-norms
1070+
of the matrices of the system pencil S.
1071+
1072+
The balancing can be performed optionally on the following
1073+
particular system pencils
1074+
1075+
S = A-lambda E,
1076+
1077+
S = ( A-lambda E B ), or
1078+
1079+
S = ( A-lambda E ).
1080+
( C )
1081+
Required arguments:
1082+
l : input int
1083+
The number of rows of matrices A, B, and E. l >= 0.
1084+
n : input int
1085+
The number of columns of matrices A, E, and C. n >= 0.
1086+
m : input int
1087+
The number of columns of matrix B. m >= 0.
1088+
p : input int
1089+
The number of rows of matrix C. P >= 0.
1090+
A : rank-2 array('d') with bounds (l,n)
1091+
The leading L-by-N part of this array must
1092+
contain the state dynamics matrix A.
1093+
E : rank-2 array('d') with bounds (l,n)
1094+
The leading L-by-N part of this array must
1095+
contain the descriptor matrix E.
1096+
B : rank-2 array('d') with bounds (l,m)
1097+
The leading L-by-M part of this array must
1098+
contain the input/state matrix B.
1099+
The array B is not referenced if M = 0.
1100+
C : rank-2 array('d') with bounds (p,n)
1101+
The leading P-by-N part of this array must
1102+
contain the state/output matrix C.
1103+
The array C is not referenced if P = 0.
1104+
Optional arguments:
1105+
job := 'A' input string(len=1)
1106+
Indicates which matrices are involved in balancing, as
1107+
follows:
1108+
= 'A': All matrices are involved in balancing;
1109+
= 'B': B, A and E matrices are involved in balancing;
1110+
= 'C': C, A and E matrices are involved in balancing;
1111+
= 'N': B and C matrices are not involved in balancing.
1112+
thresh := 0.0 input float
1113+
Threshold value for magnitude of elements:
1114+
elements with magnitude less than or equal to
1115+
THRESH are ignored for balancing. THRESH >= 0.
1116+
Return objects:
1117+
A : rank-2 array('d') with bounds (l,n)
1118+
The leading L-by-N part of this array contains
1119+
the balanced matrix Dl*A*Dr.
1120+
E : rank-2 array('d') with bounds (l,n)
1121+
The leading L-by-N part of this array contains
1122+
the balanced matrix Dl*E*Dr.
1123+
B : rank-2 array('d') with bounds (l,m)
1124+
If M > 0, the leading L-by-M part of this array
1125+
contains the balanced matrix Dl*B.
1126+
The array B is not referenced if M = 0.
1127+
C : rank-2 array('d') with bounds (p,n)
1128+
If P > 0, the leading P-by-N part of this array
1129+
contains the balanced matrix C*Dr.
1130+
The array C is not referenced if P = 0.
1131+
lscale : rank-1 array('d') with bounds (l)
1132+
The scaling factors applied to S from left. If Dl(j) is
1133+
the scaling factor applied to row j, then
1134+
SCALE(j) = Dl(j), for j = 1,...,L.
1135+
rscale : rank-1 array('d') with bounds (n)
1136+
The scaling factors applied to S from right. If Dr(j) is
1137+
the scaling factor applied to column j, then
1138+
SCALE(j) = Dr(j), for j = 1,...,N.
1139+
"""
1140+
1141+
hidden = ' (hidden by the wrapper)'
1142+
arg_list = ['job', 'l', 'n', 'm', 'p', 'thresh', 'A', 'lda'+hidden, 'E','lde'+hidden,'B','ldb'+hidden,'C','ldc'+hidden, 'lscale', 'rscale', 'dwork'+hidden, 'info']
1143+
1144+
if job != 'A' and job != 'B' and job != 'C' and job != 'N':
1145+
raise ValueError('Parameter job had an illegal value')
1146+
1147+
A,E,B,C,lscale,rscale,info = _wrapper.tg01ad(job,l,n,m,p,thresh,A,E,B,C)
1148+
1149+
if info < 0:
1150+
error_text = "The following argument had an illegal value: "+arg_list[-info-1]
1151+
e = ValueError(error_text)
1152+
e.info = info
1153+
raise e
1154+
if info != 0:
1155+
e = ArithmeticError('tg01ad failed')
1156+
e.info = info
1157+
raise e
1158+
1159+
return A,E,B,C,lscale,rscale
1160+
10531161
def tg01fd(l,n,m,p,A,E,B,C,Q=None,Z=None,compq='N',compz='N',joba='N',tol=0.0,ldwork=None):
10541162
""" A,E,B,C,ranke,rnka22,Q,Z = tg01fd(l,n,m,p,A,E,B,C,[Q,Z,compq,compz,joba,tol,ldwork])
10551163

0 commit comments

Comments
 (0)