Skip to content

Correct td04ad #27

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
murrayrm merged 2 commits into from
Jun 29, 2018
Merged

Correct td04ad #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
30dc0e5
work on getting 'C' version for td04ad working and checked
repagh Jun 4, 2018
7fce83b
working TD04AD.f, does not crash any more on edge cases (N=0, NR=0),
repagh Jun 20, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 17 additions & 6 deletions slycot/src/TD04AD.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,8 @@ SUBROUTINE TD04AD( ROWCOL, M, P, INDEX, DCOEFF, LDDCOE, UCOEFF,
$ UCOEFF(LDUCO1,LDUCO2,*)
C .. Local Scalars ..
LOGICAL LROCOC, LROCOR
INTEGER I, J, JSTOP, K, KDCOEF, MPLIM, MWORK, N, PWORK
INTEGER I, J, JSTOP, K, KDCOEF, MPLIM, MWORK, N, PWORK,
$ KU
C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
Expand Down Expand Up @@ -401,18 +402,28 @@ SUBROUTINE TD04AD( ROWCOL, M, P, INDEX, DCOEFF, LDDCOE, UCOEFF,
C state-space representation, and reorder the rows and columns
C to get an upper block Hessenberg state dynamics matrix.
C
K = IWORK(1)+IWORK(2)-1
CALL TB01XD( 'D', NR, MWORK, PWORK, K, NR-1, A, LDA, B, LDB,
$ C, LDC, D, LDD, INFO )
C IWORK contains the orders of the diagnonal blocks
C RvP, In TB01PD, IWORK is zeroed from INDCON to N, beyond N it may
C contain nonsense?
K = -1
DO 55 I = 1, N
K = K + IWORK(I)
55 CONTINUE
C
C RvP 180615 Try to protect against re-working an empty [] A
C matrix, failed with K < 0
C
CALL TB01XD( 'D', NR, MWORK, PWORK, MAX(0, K), MAX(0,NR-1),
$ A, LDA, B, LDB, C, LDC, D, LDD, INFO )
IF ( MPLIM.NE.1 ) THEN
C
C Also, retranspose U(s) if this is non-scalar.
C
DO 70 K = 1, KDCOEF
C
C
DO 60 J = 1, JSTOP
CALL DSWAP( MPLIM-J, UCOEFF(J+1,J,K), 1,
$ UCOEFF(J,J+1,K), LDUCO1 )
$ UCOEFF(J,J+1,K), LDUCO1 )
60 CONTINUE
C
70 CONTINUE
Expand Down
6 changes: 3 additions & 3 deletions slycot/src/transform.pyf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -344,17 +344,17 @@ subroutine td04ad_c(rowcol,m,p,index_bn,dcoeff,lddcoe,ucoeff,lduco1,lduco2,nr,a,
integer check(m>=0) :: m
integer check(p>=0) :: p
integer dimension(m),depend(m) :: index_bn
double precision intent(in,copy),dimension(max(1,m),*),depend(p) :: dcoeff
double precision intent(in,copy),dimension(max(1,m),*),depend(m) :: dcoeff
integer intent(hide),depend(dcoeff) :: lddcoe=shape(dcoeff,0)
double precision intent(in,copy),dimension(max(1,max(m,p)),max(1,max(m,p)),*),depend(p,m) :: ucoeff
integer intent(hide),depend(ucoeff) :: lduco1=shape(ucoeff,0)
integer intent(hide),depend(ucoeff) :: lduco2=shape(ucoeff,1)
integer intent(in,out) :: nr != sum(index_bn)
double precision intent(out),dimension(max(1,nr),nr),depend(nr) :: a
double precision intent(out),dimension(max(1,nr),max(1,nr)),depend(nr) :: a
integer intent(hide),depend(a) :: lda = shape(a,0)
double precision intent(out),dimension(max(1,nr),max(m,p)),depend(nr,m,p) :: b
integer intent(hide),depend(b) :: ldb = shape(b,0)
double precision intent(out),dimension(max(1,max(m,p)),nr),depend(nr,m,p) :: c
double precision intent(out),dimension(max(1,max(m,p)),max(1,nr)),depend(nr,m,p) :: c
integer intent(hide),depend(c) :: ldc = shape(c,0)
double precision intent(out),dimension(max(1,max(m,p)),max(m,p)),depend(p,m) :: d
integer intent(hide),depend(d) :: ldd = shape(d,0)
Expand Down
138 changes: 138 additions & 0 deletions slycot/tests/test_td04ad.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
#!/usr/bin/env python
#
# test_td04ad.py - test suite for tf -> ss conversion
# RvP, 04 Jun 2018

from __future__ import print_function

import unittest
from slycot import transform
import numpy as np

from numpy.testing import assert_raises, assert_almost_equal

class TestTf2SS(unittest.TestCase):

def test_td04ad_case1(self):
"""td04ad: Convert with both 'C' and 'R' options"""

# for octave:
"""
num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ];
[3.0, -1.0, 1.0 ], [ 0.0, 1.0 ];
[0.0, 0.0, 1.0], [ 0.0, 2.0 ] };
den = { [1.0, 0.4, 3.0], [ 1.0, 1.0 ];
[1.0, 0.4, 3.0], [ 1.0, 1.0 ];
[1.0, 0.4, 3.0], [ 1.0, 1.0 ]};
"""

# common denominators for the inputs
n = 2
m = 2
p = 3
num = np.array([
[ [0.0, 0.0, 1.0 ], [ 1.0, 0.0, 0.0 ] ],
[ [3.0, -1.0, 1.0 ], [ 0.0, 1.0, 0.0 ] ],
[ [0.0, 0.0, 1.0], [ 0.0, 2.0, 0.0 ] ] ])
p, m, d = num.shape
numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
numc[:p,:m,:] = num

denc = np.array(
[ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ])
indc = np.array(
[ 2, 1 ], dtype=int)
denr = np.array(
[ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ], [1.0, 0.0, 0.0] ])
indr = np.array(
[ 2, 1, 0 ], dtype=int)

n, A, B, C, D = transform.td04ad('C', 2, 3, indc, denc, numc)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
Ac = [ [-1, 0, 0], [ 0, -0.4, -0.3], [ 0, 10, 0]]
Bc = [ [0, -1] ,[ 1 , 0], [ 0, 0]]
Cc = [ [1, 0, 0.1], [-1, -2.2, -0.8], [ -2, 0, 0.1] ]
Dc = [ [0, 1], [ 3, 0], [ 0, 0]]
np.testing.assert_array_almost_equal(A, Ac)
np.testing.assert_array_almost_equal(B, Bc)
np.testing.assert_array_almost_equal(C, Cc)
np.testing.assert_array_almost_equal(D, Dc)

resr = transform.td04ad('R', 2, 3, indr, denr, num)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)

def test_staticgain(self):
"""td04ad: Convert a transferfunction to SS with only static gain"""

# 2 inputs, 3 outputs? columns share a denominator
num = np.array([ [ [1.0], [2.0] ],
[ [0.2], [4.3] ],
[ [1.2], [3.2] ] ])
p, m, d = num.shape
numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
numc[:p,:m,:] = num

# denc, columns share a common denominator
denc = np.array([ [ 1.0], [0.5] ])
Dc = (num / denc).reshape((3,2))
idxc = np.zeros((2,), dtype=int)

# denr, rows share a common denominator
denr = np.array([ [1.0], [0.5], [3.0] ])
idxr = np.zeros((3,), dtype=int)
Dr = (num / denr[:, np.newaxis]).reshape((3,2))

# fails with:
# On entry to TB01XD parameter number 5 had an illegal value

n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
self.assertEqual(A.shape, (0,0))
self.assertEqual(B.shape, (0,2))
self.assertEqual(C.shape, (3,0))
np.testing.assert_array_almost_equal(D, Dc)

n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num)
#print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D)
self.assertEqual(A.shape, (0,0))
self.assertEqual(B.shape, (0,2))
self.assertEqual(C.shape, (3,0))
np.testing.assert_array_almost_equal(D, Dr)


def test_mixfeedthrough(self):
"""Test case popping up from control testing"""
# a mix of feedthrough and dynamics. The problem from the control
# package was somewhere else
num = np.array([ [ [ 0.0, 0.0 ], [ 0.0, -0.2 ] ],
[ [ -0.1, 0.0 ], [ 0.0, 0.0 ] ] ])
p, m, d = num.shape
numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float)
numc[:p,:m,:] = num
denc = np.array([ [ 1.0, 1.1 ], [ 1.0, 0.0 ] ])
idxc = np.array([ 1, 0 ])
n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc)
np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]]))

def test_toandfrom(self):

A = np.array([[-3.0]])
B = np.array([[0.1, 0.0]])
C = np.array([[1.0],[0.0]])
D = np.array([[0.0, 0.0],[0.0, 1.0]])

tfout = transform.tb04ad(1, 2, 2, A, B, C, D)

num = tfout[6]
den = tfout[5]
idxc = np.array([1, 0])
n, At, Bt, Ct, Dt = transform.td04ad('R', 2, 2, idxc, den, num)
np.testing.assert_array_almost_equal(D, Dt)
np.testing.assert_array_almost_equal(A, At)

def suite():
return unittest.TestLoader().loadTestsFromTestCase(TestTF2SS)


if __name__ == "__main__":
unittest.main()
6 changes: 4 additions & 2 deletions slycot/transform.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -670,7 +670,7 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None):
e = ValueError("The numerator is not a 3D array!")
e.info = -7
raise e
if ucoeff.shape != (max([1,m,p]),max([1,m,p]),kdcoef):
if ucoeff.shape != (max(1,m,p),max(1,m,p),kdcoef):
e = ValueError("The numerator shape is ("+str(ucoeff.shape[0])+","+str(ucoeff.shape[1])+","+str(ucoeff.shape[2])+"), but expected ("+str(max([1,m,p]))+","+str(max([1,m,p]))+","+str(kdcoef)+")")
e.info = -7
raise e
Expand All @@ -690,7 +690,9 @@ def td04ad(rowcol,m,p,index,dcoeff,ucoeff,tol=0.0,ldwork=None):
e.info = out[-1]
raise e
if out[-1] > 0:
error_text = "The leading coefficient of a denominator polynomial is nearly zero; calculations would overflow; no state-space representation was calculated. ABS(DCOEFF("+str(out[-1])+",1))="+str(abs(dcoeff(out[-1],1)))+" is too small."
error_text = "The leading coefficient of a denominator polynomial is nearly zero; calculations would overflow; no state-space representation was calculated. ABS(DCOEFF("+str(out[-1])+",1))="+str(abs(dcoeff[out[-1],1]))+" is too small."
print(dcoeff)
e = ValueError(error_text)
e.info = out[-1]
raise e
Nr, A, B, C, D = out[:-1]
Expand Down