Skip to content
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.

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

Bugfix ab13bd #200

Merged
merged 4 commits into from
Aug 24, 2023
Merged

Bugfix ab13bd #200

merged 4 commits into from
Aug 24, 2023

Conversation

KybernetikJo
Copy link
Contributor

This Draft/PR fixes the issues #198 and #199.

@bnavigator
Copy link
Collaborator

Hi, thanks for your contribution. I didn't check in detail yet, but if it is really true that there is a call by reference, then this would be the case in many more functions than ab13bd and we have a general problem in the whole package.

@bnavigator bnavigator added the bug label Jul 20, 2023
@KybernetikJo
Copy link
Contributor Author

The following changes are actually enough:

    integer check(n>=0) :: n
    integer check(m>=0) :: m
    integer check(p>=0) :: p
    double precision intent(in,out,copy), dimension(n,n),depend(n) :: a
    integer intent(hide),depend(a) :: lda = shape(a,0)
    double precision intent(in,out,copy), dimension(n,m),depend(n,m) :: b
    integer intent(hide),depend(b) :: ldb = shape(b,0)
    double precision intent(in,out,copy), dimension(p,n),depend(n,p) :: c
    integer intent(hide),depend(c) :: ldc = shape(c,0)
    double precision intent(in,out,copy), dimension(p,m),depend(m,p) :: d
    integer intent(hide),depend(d) :: ldd = shape(d,0)

Intent(in,out,copy) copies the parameter. The function has no side effects on the state space matrices A,B,C,D

TODO:

  • Add unit tests, pytests
  • Further testing.
  • Doing background reading about f2py.

@bnavigator
Copy link
Collaborator

Actually, according to https://numpy.org/doc/stable/f2py/signature-file.html, we should have the desired intent(in) by default.

in
The corresponding argument is considered to be input-only. This means that the value of the argument is passed to a
Fortran/C function and that the function is expected to not change the value of this argument.

The following rules apply:

  • If none of intent(in | inout | out | hide) are specified, intent(in) is assumed.

I wonder if it the stray

double precision intent(out) :: ab13bd

causes the inputs to be altered?

Could you check if you get altered inputs in other wrappers as well?

@bnavigator
Copy link
Collaborator

Intent(in,out,copy) copies the parameter.

Yes, but we actually don't want the output anywhere.

@KybernetikJo
Copy link
Contributor Author

Thanks for your reply. I appreciate that.

First, I am going to build up some theoretical understanding about f2py. Then I can run some experiments.

Intent(in,out,copy) copies the parameter.

Yes, but we actually don't want the output anywhere.

I am not sure. Why SLICOT returns them?

@KybernetikJo KybernetikJo mentioned this pull request Jul 20, 2023
3 tasks
@roryyorke
Copy link
Collaborator

FWIW I can reproduce this with Slycot 0.5.2 on Ubuntu 20.04 in a Conda environment.

From the SLICOT docs for AB13MD, the A, B, C and D arguments are input/output; they contain the state-space matrices of "the numerator factor Q of the right coprime factorization with inner denominator of G". So I think we violate this expectation of f2py:

in
The corresponding argument is considered to be input-only. This means that the value of the argument is passed to a
Fortran/C function and that the function is expected to not change the value of this argument.

Specifically, we violate "the function is expected to not change the value of this argument". Whoops! As @bnavigator noted, this could be a problem in quite a few places, we'll have to check.

@roryyorke
Copy link
Collaborator

See https://github.com/roryyorke/Slycot/tree/ab13md-alter-args, whence:

  • 0277e72 - add regression test from @KybernetikJo 's example; test fails

  • 4699356 , remove (strange?) intent(out) :: ab13bd, test fails:

diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
index 633ff6ec..183a22f2 100644
--- a/slycot/src/analysis.pyf
+++ b/slycot/src/analysis.pyf
@@ -321,7 +321,6 @@ function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwar
     integer intent(hide),depend(n,m,p) :: ldwork = max(1,max(m*(n+m)+max(n*(n+5),max(m*(m+2),4*p)),n*(max(n,p)+4)+min(n,p)))
     integer intent(out) :: iwarn
     integer intent(out) :: info
-    double precision intent(out) :: ab13bd
 end function ab13bd
 subroutine ab13dd(dico,jobe,equil,jobd,n,m,p,fpeak,a,lda,e,lde,b,ldb,c,ldc,d,ldd,gpeak,tol,iwork,dwork,ldwork,cwork,lcwork,info) ! in AB13DD.f
     character intent(in) :: dico
diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
index 183a22f2..999e3b65 100644
--- a/slycot/src/analysis.pyf
+++ b/slycot/src/analysis.pyf
@@ -307,13 +307,13 @@ function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwar
     integer check(n>=0) :: n
     integer check(n>=0) :: m
     integer check(n>=0) :: p
-    double precision dimension(n,n),depend(n) :: a
+    double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
     integer intent(hide),depend(a) :: lda = shape(a,0)
-    double precision dimension(n,m),depend(n,m) :: b
+    double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b
     integer intent(hide),depend(b) :: ldb = shape(b,0)
-    double precision dimension(p,n),depend(n,p) :: c
+    double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c
     integer intent(hide),depend(c) :: ldc = shape(c,0)
-    double precision dimension(p,m),depend(m,p) :: d
+    double precision intent(in,out,copy),dimension(p,m),depend(m,p) :: d
     integer intent(hide),depend(d) :: ldd = shape(d,0)
     integer intent(out) :: nq
     double precision :: tol

@roryyorke
Copy link
Collaborator

and 2625e1f, intent(in,copy) also passes.

@roryyorke
Copy link
Collaborator

I don't find the f2py docs all that clear. I don't understand the difference between intent variants in,out, inout, and inplace; nor the difference between copy and overwrite.

For the record, intent(in, overwrite) fails the ab13bd test.

@roryyorke
Copy link
Collaborator

according to

grep -n "dimension"  ~/src/slycot/slycot/src/*.pyf|grep -v "intent" >> audit.out
grep -n "dimension" ~/src/slycot/slycot/src/*.pyf|fgrep "intent(in)"  >> audit.out

there are 139 array arguments with implicit or explicit intent(in) on master (including the 4 identified by @KybernetikJo ).

grep results
/home/rory/src/slycot/slycot/src/analysis.pyf:30:    double precision dimension(n1,n1),depend(n1) :: a1
/home/rory/src/slycot/slycot/src/analysis.pyf:32:    double precision dimension(n1,m1),depend(n1,m1) :: b1
/home/rory/src/slycot/slycot/src/analysis.pyf:34:    double precision dimension(p1,n1),depend(n1,p1) :: c1
/home/rory/src/slycot/slycot/src/analysis.pyf:36:    double precision dimension(p1,m1),depend(m1,p1) :: d1
/home/rory/src/slycot/slycot/src/analysis.pyf:38:    double precision dimension(n2,n2),depend(n2) :: a2
/home/rory/src/slycot/slycot/src/analysis.pyf:40:    double precision dimension(n2,p1),depend(n2,p1) :: b2
/home/rory/src/slycot/slycot/src/analysis.pyf:42:    double precision dimension(p2,n2),depend(n2,p2) :: c2
/home/rory/src/slycot/slycot/src/analysis.pyf:44:    double precision dimension(p2,p1),depend(p1,p2) :: d2
/home/rory/src/slycot/slycot/src/analysis.pyf:66:    double precision dimension(n1,n1),depend(n1) :: a1
/home/rory/src/slycot/slycot/src/analysis.pyf:68:    double precision dimension(n1,m1),depend(n1,m1) :: b1
/home/rory/src/slycot/slycot/src/analysis.pyf:70:    double precision dimension(p1,n1),depend(p1,n1) :: c1
/home/rory/src/slycot/slycot/src/analysis.pyf:72:    double precision dimension(p1,m1),depend(p1,m1) :: d1
/home/rory/src/slycot/slycot/src/analysis.pyf:74:    double precision dimension(n2,n2),depend(n2) :: a2
/home/rory/src/slycot/slycot/src/analysis.pyf:76:    double precision dimension(n2,p1),depend(n2,p1) :: b2
/home/rory/src/slycot/slycot/src/analysis.pyf:78:    double precision dimension(m1,n2),depend(m1,n2) :: c2
/home/rory/src/slycot/slycot/src/analysis.pyf:80:    double precision dimension(m1,p1),depend(m1,p1) :: d2
/home/rory/src/slycot/slycot/src/analysis.pyf:310:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/analysis.pyf:312:    double precision dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/analysis.pyf:314:    double precision dimension(p,n),depend(n,p) :: c
/home/rory/src/slycot/slycot/src/analysis.pyf:316:    double precision dimension(p,m),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/analysis.pyf:335:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/analysis.pyf:337:    double precision dimension(n,n),depend(n) :: e
/home/rory/src/slycot/slycot/src/analysis.pyf:339:    double precision dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/analysis.pyf:341:    double precision dimension(p,n),depend(n,p) :: c
/home/rory/src/slycot/slycot/src/analysis.pyf:343:    double precision dimension(p,m),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/analysis.pyf:375:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/analysis.pyf:386:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:10:    double precision dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:37:	double precision dimension(n,n),depend(n) :: g
/home/rory/src/slycot/slycot/src/synthesis.pyf:68:	double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:70:	double precision dimension(n,m),depend(m,n) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:72:	double precision dimension(n,n),depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:74:	double precision dimension(m,m),depend(m) :: r
/home/rory/src/slycot/slycot/src/synthesis.pyf:76:	double precision dimension(n,m),depend(m,n) :: l
/home/rory/src/slycot/slycot/src/synthesis.pyf:108:	double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:110:	double precision dimension(n,m),depend(m,n) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:112:	double precision dimension(p,n),depend(n,p) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:114:	double precision dimension(m,m),depend(m) :: r
/home/rory/src/slycot/slycot/src/synthesis.pyf:116:	double precision dimension(n,m),depend(m,n) :: l
/home/rory/src/slycot/slycot/src/synthesis.pyf:148:	double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:150:	double precision dimension(n,m),depend(m,n) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:152:	double precision dimension(n,n),depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:154:	double precision dimension(p,m),depend(m,p) :: r
/home/rory/src/slycot/slycot/src/synthesis.pyf:156:	double precision dimension(n,m),depend(m,n) :: l
/home/rory/src/slycot/slycot/src/synthesis.pyf:188:	double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:190:	double precision dimension(n,m),depend(m,n) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:192:	double precision dimension(p,n),depend(n,p) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:194:	double precision dimension(p,m),depend(m,p) :: r
/home/rory/src/slycot/slycot/src/synthesis.pyf:196:	double precision dimension(n,m),depend(m,n) :: l
/home/rory/src/slycot/slycot/src/synthesis.pyf:358:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:360:    double precision dimension(n,n),depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:417:    double precision dimension(n,n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:419:    double precision dimension(n,m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:421:    double precision dimension(np,n), depend(np,n) :: c
/home/rory/src/slycot/slycot/src/synthesis.pyf:423:    double precision dimension(np,m), depend(np,m) :: d
/home/rory/src/slycot/slycot/src/synthesis.pyf:464:    double precision dimension(n,n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:466:    double precision dimension(n,m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:468:    double precision dimension(np,n), depend(np,n) :: c
/home/rory/src/slycot/slycot/src/synthesis.pyf:470:    double precision dimension(np,m), depend(np,m) :: d
/home/rory/src/slycot/slycot/src/synthesis.pyf:535:    double precision dimension(n,n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:537:    double precision dimension(n,m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:539:    double precision dimension(np,n), depend(np,n) :: c
/home/rory/src/slycot/slycot/src/synthesis.pyf:541:    double precision dimension(np,m), depend(np,m) :: d
/home/rory/src/slycot/slycot/src/synthesis.pyf:619:    double precision dimension(max(1,n),n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:621:    double precision dimension(max(1,n),n), depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:623:    double precision dimension(max(1,n),n), depend(n) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:625:    double precision dimension(max(1,n),n), depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:664:    double precision dimension(max(1,n),n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:666:    double precision dimension(max(1,n),n), depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:668:    double precision dimension(max(1,n),m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:670:    double precision dimension(max(1,n),n), depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:672:    double precision dimension(max(1,m),m), depend(m) :: r ! fact = 'N'
/home/rory/src/slycot/slycot/src/synthesis.pyf:674:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
/home/rory/src/slycot/slycot/src/synthesis.pyf:709:    double precision dimension(max(1,n),n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:711:    double precision dimension(max(1,n),n), depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:713:    double precision dimension(max(1,n),m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:715:    double precision dimension(max(1,p),n), depend(p,n) :: q ! fact = 'C'
/home/rory/src/slycot/slycot/src/synthesis.pyf:717:    double precision dimension(max(1,m),m), depend(m) :: r ! fact = 'C'
/home/rory/src/slycot/slycot/src/synthesis.pyf:719:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
/home/rory/src/slycot/slycot/src/synthesis.pyf:754:    double precision dimension(max(1,n),n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:756:    double precision dimension(max(1,n),n), depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:758:    double precision dimension(max(1,n),m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:760:    double precision dimension(max(1,n),n), depend(n) :: q ! fact = 'D'
/home/rory/src/slycot/slycot/src/synthesis.pyf:762:    double precision dimension(max(1,p),m), depend(p,m) :: r ! fact = 'D'
/home/rory/src/slycot/slycot/src/synthesis.pyf:764:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
/home/rory/src/slycot/slycot/src/synthesis.pyf:799:    double precision dimension(max(1,n),n), depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:801:    double precision dimension(max(1,n),n), depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:803:    double precision dimension(max(1,n),m), depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:805:    double precision dimension(max(1,p),n), depend(p,n) :: q ! fact = 'B'
/home/rory/src/slycot/slycot/src/synthesis.pyf:807:    double precision dimension(max(1,p),m), depend(p,m) :: r ! fact = 'B'
/home/rory/src/slycot/slycot/src/synthesis.pyf:809:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
/home/rory/src/slycot/slycot/src/synthesis.pyf:838:    double precision dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:840:    double precision dimension(n,n),depend(n) :: e
/home/rory/src/slycot/slycot/src/synthesis.pyf:842:    double precision dimension(n,n),depend(n) :: q
/home/rory/src/slycot/slycot/src/synthesis.pyf:844:    double precision dimension(n,n),depend(n) :: z
/home/rory/src/slycot/slycot/src/transform.pyf:32:	double precision dimension(max(m,p),max(m,p)),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:64:	double precision dimension(max(m,p),max(m,p)),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:253:	double precision dimension(m,m,indlim),depend(p,indlim) :: pcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:256:	double precision dimension(max(m,p),max(m,p),indlim),depend(m,p,indlim) :: qcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:266:	integer dimension(p) :: index_bn
/home/rory/src/slycot/slycot/src/transform.pyf:267:	double precision required,dimension(p,p,*),depend(p) :: pcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:270:	double precision required,dimension(p,m,*),depend(m,p) :: qcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:293:	integer dimension(m) :: index_bn
/home/rory/src/slycot/slycot/src/transform.pyf:294:	double precision required,dimension(m,m,*),depend(m) :: pcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:297:	double precision required,dimension(max(m,p),max(m,p),*),depend(m,p) :: qcoeff
/home/rory/src/slycot/slycot/src/transform.pyf:320:    integer dimension(p),depend(p) :: index_bn
/home/rory/src/slycot/slycot/src/transform.pyf:346:    integer dimension(m),depend(m) :: index_bn
/home/rory/src/slycot/slycot/src/transform.pyf:372:    double precision required,dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/transform.pyf:374:    double precision required,dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/transform.pyf:376:    double precision required,dimension(p,n),depend(n,p) :: c
/home/rory/src/slycot/slycot/src/transform.pyf:378:    double precision required,dimension(p,m),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:380:    double precision required,dimension(m,ny),depend(m,ny) :: u
/home/rory/src/slycot/slycot/src/transform.pyf:393:    double precision dimension(na,na),depend(na) :: a
/home/rory/src/slycot/slycot/src/transform.pyf:395:    double precision dimension(na,nb),depend(na,nb) :: b
/home/rory/src/slycot/slycot/src/transform.pyf:397:    double precision dimension(nc,na),depend(nc,na) :: c
/home/rory/src/slycot/slycot/src/analysis.pyf:118:    double precision intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a
/home/rory/src/slycot/slycot/src/analysis.pyf:120:    double precision intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b
/home/rory/src/slycot/slycot/src/analysis.pyf:122:    double precision intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c
/home/rory/src/slycot/slycot/src/analysis.pyf:124:    double precision intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d
/home/rory/src/slycot/slycot/src/analysis.pyf:149:    complex*16 intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a
/home/rory/src/slycot/slycot/src/analysis.pyf:151:    complex*16 intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b
/home/rory/src/slycot/slycot/src/analysis.pyf:153:    complex*16 intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c
/home/rory/src/slycot/slycot/src/analysis.pyf:155:    complex*16 intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d
/home/rory/src/slycot/slycot/src/analysis.pyf:357:  complex*16 intent(in),dimension(n,n) :: z
/home/rory/src/slycot/slycot/src/analysis.pyf:360:  integer intent(in),dimension(m) :: nblock
/home/rory/src/slycot/slycot/src/analysis.pyf:361:  integer intent(in),dimension(m), depend(m) :: itype
/home/rory/src/slycot/slycot/src/math.pyf:7:	double precision intent(in),check(shape(p,0)==dp+1),dimension(dp+1),depend(dp) :: p
/home/rory/src/slycot/slycot/src/math.pyf:55:    double precision intent(in),depend(n,p),dimension(ldtau,p) :: tau
/home/rory/src/slycot/slycot/src/math.pyf:105:  double precision intent(in),dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:500:    double precision intent(in),dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/synthesis.pyf:502:    double precision intent(in),dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/synthesis.pyf:504:    double precision intent(in),dimension(np,n),depend(np,n) :: c
/home/rory/src/slycot/slycot/src/synthesis.pyf:506:    double precision intent(in),dimension(np,m),depend(np,m) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:95:    double precision intent(in),dimension(max(1,p),m),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:123:    double precision intent(in),dimension(max(1,max(m,p)),max(m,p)),depend(m,p) :: d
/home/rory/src/slycot/slycot/src/transform.pyf:211:  double precision intent(in),dimension(n,n),depend(n) :: a
/home/rory/src/slycot/slycot/src/transform.pyf:213:  double precision intent(in),dimension(n,m),depend(n,m) :: b
/home/rory/src/slycot/slycot/src/transform.pyf:215:  double precision intent(in),dimension(p,n),depend(n,p) :: c

Some of these are fine, e.g., the first arg to AB05MD is marked SLICOT as "input".

We could wholesale change all of these to intent(in,copy), but it's probably best to audit all these.

Worst-case is that some dependency is relying on this bad behaviour, but I think that is unlikely.

ab13bd specifically is not used in python-control, at least.

@roryyorke
Copy link
Collaborator

I've compared each case from above with the comments in the associated SLICOT files, and marked each line OK or NOTOK. There are 19 "NOTOK" lines, including the 4 already found.

I'd appreciate spot checks of the results.

Next steps:

  • Identify Slycot problem functions. I think we are ok in any case which:
    • SLICOT args are internal to Slycot routines (I doubt there are many of these)
    • if any Slycot args are regarded as input-output, it's possibly also not a bug (as in, does not produce undesirable behaviour)
    • the idea is to prioritise user-visible cases, but eventually fix all cases
  • which Slycot problem functions are used in python-control?
SLICOT-checked use of "intent(in)" arguments
OK /home/rory/src/slycot/slycot/src/analysis.pyf:30:    double precision dimension(n1,n1),depend(n1) :: a1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:32:    double precision dimension(n1,m1),depend(n1,m1) :: b1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:34:    double precision dimension(p1,n1),depend(n1,p1) :: c1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:36:    double precision dimension(p1,m1),depend(m1,p1) :: d1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:38:    double precision dimension(n2,n2),depend(n2) :: a2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:40:    double precision dimension(n2,p1),depend(n2,p1) :: b2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:42:    double precision dimension(p2,n2),depend(n2,p2) :: c2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:44:    double precision dimension(p2,p1),depend(p1,p2) :: d2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:66:    double precision dimension(n1,n1),depend(n1) :: a1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:68:    double precision dimension(n1,m1),depend(n1,m1) :: b1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:70:    double precision dimension(p1,n1),depend(p1,n1) :: c1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:72:    double precision dimension(p1,m1),depend(p1,m1) :: d1
OK /home/rory/src/slycot/slycot/src/analysis.pyf:74:    double precision dimension(n2,n2),depend(n2) :: a2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:76:    double precision dimension(n2,p1),depend(n2,p1) :: b2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:78:    double precision dimension(m1,n2),depend(m1,n2) :: c2
OK /home/rory/src/slycot/slycot/src/analysis.pyf:80:    double precision dimension(m1,p1),depend(m1,p1) :: d2
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:310:    double precision dimension(n,n),depend(n) :: a
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:312:    double precision dimension(n,m),depend(n,m) :: b
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:314:    double precision dimension(p,n),depend(n,p) :: c
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:316:    double precision dimension(p,m),depend(m,p) :: d
OK /home/rory/src/slycot/slycot/src/analysis.pyf:335:    double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/analysis.pyf:337:    double precision dimension(n,n),depend(n) :: e
OK /home/rory/src/slycot/slycot/src/analysis.pyf:339:    double precision dimension(n,m),depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/analysis.pyf:341:    double precision dimension(p,n),depend(n,p) :: c
OK /home/rory/src/slycot/slycot/src/analysis.pyf:343:    double precision dimension(p,m),depend(m,p) :: d
OK /home/rory/src/slycot/slycot/src/analysis.pyf:375:    double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/analysis.pyf:386:    double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:10:    double precision dimension(n,m),depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:37:	double precision dimension(n,n),depend(n) :: g
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:68:	double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:70:	double precision dimension(n,m),depend(m,n) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:72:	double precision dimension(n,n),depend(n) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:74:	double precision dimension(m,m),depend(m) :: r
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:76:	double precision dimension(n,m),depend(m,n) :: l
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:108:	double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:110:	double precision dimension(n,m),depend(m,n) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:112:	double precision dimension(p,n),depend(n,p) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:114:	double precision dimension(m,m),depend(m) :: r
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:116:	double precision dimension(n,m),depend(m,n) :: l
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:148:	double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:150:	double precision dimension(n,m),depend(m,n) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:152:	double precision dimension(n,n),depend(n) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:154:	double precision dimension(p,m),depend(m,p) :: r
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:156:	double precision dimension(n,m),depend(m,n) :: l
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:188:	double precision dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:190:	double precision dimension(n,m),depend(m,n) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:192:	double precision dimension(p,n),depend(n,p) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:194:	double precision dimension(p,m),depend(m,p) :: r
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:196:	double precision dimension(n,m),depend(m,n) :: l
NOTOK /home/rory/src/slycot/slycot/src/synthesis.pyf:358:    double precision dimension(n,n),depend(n) :: a
NOTOK? /home/rory/src/slycot/slycot/src/synthesis.pyf:360:    double precision dimension(n,n),depend(n) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:417:    double precision dimension(n,n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:419:    double precision dimension(n,m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:421:    double precision dimension(np,n), depend(np,n) :: c
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:423:    double precision dimension(np,m), depend(np,m) :: d
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:464:    double precision dimension(n,n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:466:    double precision dimension(n,m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:468:    double precision dimension(np,n), depend(np,n) :: c
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:470:    double precision dimension(np,m), depend(np,m) :: d
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:535:    double precision dimension(n,n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:537:    double precision dimension(n,m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:539:    double precision dimension(np,n), depend(np,n) :: c
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:541:    double precision dimension(np,m), depend(np,m) :: d
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:619:    double precision dimension(max(1,n),n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:621:    double precision dimension(max(1,n),n), depend(n) :: e
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:623:    double precision dimension(max(1,n),n), depend(n) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:625:    double precision dimension(max(1,n),n), depend(n) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:664:    double precision dimension(max(1,n),n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:666:    double precision dimension(max(1,n),n), depend(n) :: e
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:668:    double precision dimension(max(1,n),m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:670:    double precision dimension(max(1,n),n), depend(n) :: q
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:672:    double precision dimension(max(1,m),m), depend(m) :: r ! fact = 'N'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:674:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:709:    double precision dimension(max(1,n),n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:711:    double precision dimension(max(1,n),n), depend(n) :: e
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:713:    double precision dimension(max(1,n),m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:715:    double precision dimension(max(1,p),n), depend(p,n) :: q ! fact = 'C'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:717:    double precision dimension(max(1,m),m), depend(m) :: r ! fact = 'C'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:719:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:754:    double precision dimension(max(1,n),n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:756:    double precision dimension(max(1,n),n), depend(n) :: e
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:758:    double precision dimension(max(1,n),m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:760:    double precision dimension(max(1,n),n), depend(n) :: q ! fact = 'D'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:762:    double precision dimension(max(1,p),m), depend(p,m) :: r ! fact = 'D'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:764:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:799:    double precision dimension(max(1,n),n), depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:801:    double precision dimension(max(1,n),n), depend(n) :: e
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:803:    double precision dimension(max(1,n),m), depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:805:    double precision dimension(max(1,p),n), depend(p,n) :: q ! fact = 'B'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:807:    double precision dimension(max(1,p),m), depend(p,m) :: r ! fact = 'B'
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:809:    double precision dimension(max(1,n),m), depend(n,m) :: l  ! It would be nicer to not force the dimension here when not used
NOTOK? check FACT /home/rory/src/slycot/slycot/src/synthesis.pyf:838:    double precision dimension(n,n),depend(n) :: a
NOTOK? check FACT /home/rory/src/slycot/slycot/src/synthesis.pyf:840:    double precision dimension(n,n),depend(n) :: e
NOTOK? check FACT /home/rory/src/slycot/slycot/src/synthesis.pyf:842:    double precision dimension(n,n),depend(n) :: q
NOTOK? check FACT /home/rory/src/slycot/slycot/src/synthesis.pyf:844:    double precision dimension(n,n),depend(n) :: z
OK /home/rory/src/slycot/slycot/src/transform.pyf:32:	double precision dimension(max(m,p),max(m,p)),depend(m,p) :: d
OK /home/rory/src/slycot/slycot/src/transform.pyf:64:	double precision dimension(max(m,p),max(m,p)),depend(m,p) :: d
NOTOK /home/rory/src/slycot/slycot/src/transform.pyf:253:	double precision dimension(m,m,indlim),depend(p,indlim) :: pcoeff
NOTOK /home/rory/src/slycot/slycot/src/transform.pyf:256:	double precision dimension(max(m,p),max(m,p),indlim),depend(m,p,indlim) :: qcoeff
OK /home/rory/src/slycot/slycot/src/transform.pyf:266:	integer dimension(p) :: index_bn
OK OK /home/rory/src/slycot/slycot/src/transform.pyf:267:	double precision required,dimension(p,p,*),depend(p) :: pcoeff
OK /home/rory/src/slycot/slycot/src/transform.pyf:270:	double precision required,dimension(p,m,*),depend(m,p) :: qcoeff
OK /home/rory/src/slycot/slycot/src/transform.pyf:293:	integer dimension(m) :: index_bn
OK /home/rory/src/slycot/slycot/src/transform.pyf:294:	double precision required,dimension(m,m,*),depend(m) :: pcoeff
OK /home/rory/src/slycot/slycot/src/transform.pyf:297:	double precision required,dimension(max(m,p),max(m,p),*),depend(m,p) :: qcoeff
OK /home/rory/src/slycot/slycot/src/transform.pyf:320:    integer dimension(p),depend(p) :: index_bn
OK /home/rory/src/slycot/slycot/src/transform.pyf:346:    integer dimension(m),depend(m) :: index_bn
OK /home/rory/src/slycot/slycot/src/transform.pyf:372:    double precision required,dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/transform.pyf:374:    double precision required,dimension(n,m),depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/transform.pyf:376:    double precision required,dimension(p,n),depend(n,p) :: c
OK /home/rory/src/slycot/slycot/src/transform.pyf:378:    double precision required,dimension(p,m),depend(m,p) :: d
OK /home/rory/src/slycot/slycot/src/transform.pyf:380:    double precision required,dimension(m,ny),depend(m,ny) :: u
OK /home/rory/src/slycot/slycot/src/transform.pyf:393:    double precision dimension(na,na),depend(na) :: a
OK /home/rory/src/slycot/slycot/src/transform.pyf:395:    double precision dimension(na,nb),depend(na,nb) :: b
OK /home/rory/src/slycot/slycot/src/transform.pyf:397:    double precision dimension(nc,na),depend(nc,na) :: c
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:118:    double precision intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:120:    double precision intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:122:    double precision intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c
NOTOK /home/rory/src/slycot/slycot/src/analysis.pyf:124:    double precision intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d
OK /home/rory/src/slycot/slycot/src/analysis.pyf:149:    complex*16 intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a
OK /home/rory/src/slycot/slycot/src/analysis.pyf:151:    complex*16 intent(in),dimension(ldb,*),check(shape(b,1)>=m) :: b
OK /home/rory/src/slycot/slycot/src/analysis.pyf:153:    complex*16 intent(in),dimension(ldc,*),check(shape(c,1)>=n) :: c
OK /home/rory/src/slycot/slycot/src/analysis.pyf:155:    complex*16 intent(in),dimension(ldd,*),check(shape(d,1)>=m) :: d
OK /home/rory/src/slycot/slycot/src/analysis.pyf:357:  complex*16 intent(in),dimension(n,n) :: z
OK /home/rory/src/slycot/slycot/src/analysis.pyf:360:  integer intent(in),dimension(m) :: nblock
OK /home/rory/src/slycot/slycot/src/analysis.pyf:361:  integer intent(in),dimension(m), depend(m) :: itype
OK /home/rory/src/slycot/slycot/src/math.pyf:7:	double precision intent(in),check(shape(p,0)==dp+1),dimension(dp+1),depend(dp) :: p
OK /home/rory/src/slycot/slycot/src/math.pyf:55:    double precision intent(in),depend(n,p),dimension(ldtau,p) :: tau
OK /home/rory/src/slycot/slycot/src/math.pyf:105:  double precision intent(in),dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:500:    double precision intent(in),dimension(n,n),depend(n) :: a
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:502:    double precision intent(in),dimension(n,m),depend(n,m) :: b
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:504:    double precision intent(in),dimension(np,n),depend(np,n) :: c
OK /home/rory/src/slycot/slycot/src/synthesis.pyf:506:    double precision intent(in),dimension(np,m),depend(np,m) :: d
OK /home/rory/src/slycot/slycot/src/transform.pyf:95:    double precision intent(in),dimension(max(1,p),m),depend(m,p) :: d
OK /home/rory/src/slycot/slycot/src/transform.pyf:123:    double precision intent(in),dimension(max(1,max(m,p)),max(m,p)),depend(m,p) :: d
NOTOK /home/rory/src/slycot/slycot/src/transform.pyf:211:  double precision intent(in),dimension(n,n),depend(n) :: a
NOTOK /home/rory/src/slycot/slycot/src/transform.pyf:213:  double precision intent(in),dimension(n,m),depend(n,m) :: b
NOTOK /home/rory/src/slycot/slycot/src/transform.pyf:215:  double precision intent(in),dimension(p,n),depend(n,p) :: c

@roryyorke
Copy link
Collaborator

I think there no problems at python-control level.

Fix options:

  • API-preserving and easiest: change to intent(in,out) and update Slycot docstrings where needed.
  • API-changing 1: change to intent(in,out,copy), and returned the resulting matrices
  • API-changing 2: change to intent(in, copy): users expecting mutation will be surprised.
    • removes capability for users wanting SLICOT results

Problem functions

Summary

Function Slycot args altered Python-control
ab13bd Yes Not used
sb03od Yes Used, but args local to func
sg03od Yes Not used
tc01od_r Yes Not used.
ab08nd No - mistake in analysis Used, but no problem
tb05ad_nh No - case NG alters in SLICOT, is properly handled Used, but no problem

ab13bd

H2 / L2 norm

ABCD passed straight through

docstring does not indicate that arg will be changed

not used in python-control

Conclusions:

  • needs Slycot fix
  • no python-control problem

sb03od

lyapunov equation

problem SLICOT args: a, q if FACT is not 'F'. FACT can be specified by
caller, is default 'N'.

Advertised as "in-out" in doc string.

used in python-control statefbk.gram. Args A and Q are constructed
inside the function.

Conclusions:

  • needs Slycot fix
  • no python-control problem

sg03bd

cholesky factor of generalized lyapunov equation solution

problem args: a, e, q, z. Passed straight through to SLICOT.

Not used in python-control

Conclusions:

  • needs Slycot fix
  • no python-control problem

tc01od_r

problem args: pcoeff, qcoeff

This is in,out for tc01od_l

pcoeff and qcoeff are both mutated and returned.

Conclusions

  • needs fix; not sure what
  • no python-control problem

ab08nd

mistake in initial analysis; args are input only in SLICOT

tb05ad_nh

used in slycot.tb05ad

problem args a,b,c

Complex freq. resp for A hessenberg (job='NH') A,B,C passed straight
through.

Ah - but only output if inita is "NG". For this case tb05ad_ng, with
`intent(in,out,copy)` is used.

Conclusions:

  • no bug: for NH: A,B,C not changed

@KybernetikJo
Copy link
Contributor Author

KybernetikJo commented Jul 23, 2023

I think the line

double precision intent(out) :: ab13bd
4699356

is there because of ab13bd is a function, and not a subroutine. Fortran functions do have a return value, subrountines do not.

I think ab13bd is the only function in SLICOT, I do not know why, but there might be a reason.
I am not sure deleting the line, is the right way. It might be it though.

@KybernetikJo
Copy link
Contributor Author

KybernetikJo commented Jul 23, 2023

I was wrong. ab13cd is another function. There are a few of them in SLICOT.

It also has a return value.

AB13CD = ( GAMMAL + GAMMAU )/TWO

@KybernetikJo
Copy link
Contributor Author

FWIW, I have done some background reading about F2PY, and this is my current understanding. I am a newbie to fortran, therefore you should take this with a bit of toast.

The stable version of f2py stable do have some duplication in the documentation.
I fixed that, so the dev version f2py dev is a better read.
The online "book" python fortran can also be helpful.
The f2py was design with Fortran90/95 in mind.
The goal of F2PY is, to make the use of Fortran in Python easy, and to make the API PYTHONIC. The behavior of F2PY depends heavily on the data types used in the Fortran procedures, the data types used in Python and on the F2PY signature (none, infile.f or .pyf). Therefore, very simple conclusions cannot be drawn, some understanding of Fortran seems to be necessary and good templates might be helpful.

First, Fortran has two procedures, subroutines and functions. (This effects the F2PY signature file pyf)
Both subroutines and functions have access to variables in the parent scope by argument association,
this is similar to call by reference.

All arguments in Fortran77 can be seen as call by reference, by default.

In Fortran90/95 the intent attribute has been introduced intent(in), intent(out), intent(inout). This allows the compiler to check for unintentional errors and provides self-documentation (in Fortran77 the API was only documented in comments). In addition, some compilers can do code optimization based on intent, mabye?. Is there performance advantage of using intent(in) and intent(out)?

intent in Fortran90/95
Writing fast Fortran routines for Python

intent(in): The variable is an input to the subroutine only. Its value
must not be changed during the course of the subroutine.

intent(out): The variable is an output from the subroutine only. Its
input value is irrelevant. We must assign this variable a
value before exiting the subroutine.

intent(inout): The subroutine both uses and modifies the data in the
variable. Thus, the initial value is sent and we ultimately
make modifications base on it.

In Fortran there is also the concept of a pure function, which do not have side effects. What is a pure function?
In Fortran90/95 a new keyword pure was introduced.

In Fortran2003 a value attribute has been introduced. Now call by value is possible. F2PY do not have a good support of Fortran2003 features, so it can be mostly ignored.

@KybernetikJo
Copy link
Contributor Author

KybernetikJo commented Jul 27, 2023

There are many option how to design F2PY signature files. The intent attribute in signature files as tons of options.

===================

However, currently we need a solution for the case, where the parameter are labeled as input/output in SLICOT,
are arrays in fortran and numpy arrays in python. Mostly for vectors and matrices e.g state space matrices A,B,C,D.

#200 (comment)

Fix options:

  • API-preserving and easiest: change to intent(in,out) and update Slycot docstrings where needed.

  • API-changing 1: change to intent(in,out,copy), and returned the resulting matrices

  • API-changing 2: change to intent(in, copy): users expecting mutation will be surprised.

    • removes capability for users wanting SLICOT results
  • API-changing 2: change to intent(in, copy): users expecting mutation will be surprised.

This option would potentially severely limit SLICOT procedures.
I think, if we have an INPUT/OUTPUT argument in SLICOT, we should use intent(in,out) or intent(in,out,copy).

  • API-preserving and easiest: change to intent(in,out) and update Slycot docstrings where needed.

Here, we would have to deal with Slycot function side effects somewhere.
We could use shallow copies in the outer wrapper python files.
Or we could deal with the side effects in python-control.
This could lead to many bugs.

Having said that, there are inplace function in python

[].sort()

I think in numpy there are not that common. For instance

a_out = np.sort(a)
  • API-changing 1: change to intent(in,out,copy), and returned the resulting matrices

I still think this is the best option for us. This makes the slycot._wrapper function already more PYTHONIC.

===================

We could use even more intent attribute in the signature files to make the slycot._wrapper functions easier and save to use, and more PYTHONIC.
This would make massive changes in pyhton-control necessary.

The output

double precision intent(out) :: ab13bd

in the .pyf file is there because ab13bd is function and not a subrountine. (see #200 (comment))

In the following example the optional attribute with default values is used heavily. The dimension n,m,p are determind by shape.

function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwarn,info) ! in AB13BD.f
    character optional :: dico = 'C'
    character optional :: jobn = 'H'
    integer optional :: n = shape(a,0)
    integer optional :: m = shape(b,1)
    integer optional :: p = shape(c,0)
    double precision dimension(n,n), intent(in,out,copy) :: a
    integer optional, depend(a) :: lda = shape(a,0)
    double precision dimension(n,m), intent(in,out,copy) :: b
    integer optional, depend(b) :: ldb = shape(b,0)
    double precision dimension(p,n), intent(in,out,copy) :: c
    integer optional, depend(c) :: ldc = shape(c,0)
    double precision dimension(p,m), intent(in,out,copy) :: d
    integer optional, depend(d) :: ldd = shape(d,0)
    integer intent(out) :: nq
    double precision optional :: tol = 0.0
    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
    integer optional :: ldwork = max(1,max(m*(n+m)+max(n*(n+5),max(m*(m+2),4*p)),n*(max(n,p)+4)+min(n,p)))
    integer intent(out) :: iwarn
    integer intent(out) :: info
    double precision intent(out) :: ab13bd
end function ab13bd

This would allow a minimal function call

slycot._wrapper.ab13bd(A,B,C,D)

but also some optional attributes

out = slycot._wrapper.ab13bd(A, B, C, D, dico=dico, jobn=jobn, n=n, m=m, p=p, tol=tol)

are possible.

@KybernetikJo
Copy link
Contributor Author

KybernetikJo commented Jul 27, 2023

To sum it up,

I would keep these commits and I would go for intent(in,out,copy) for the short term.

In the long run, we might reconsider the whole slycot API,
since the inner and outer wrappers for different procedures are quite heterogeneous.

Good templates could be a good solution.

@KybernetikJo KybernetikJo marked this pull request as ready for review July 27, 2023 20:42
@KybernetikJo
Copy link
Contributor Author

KybernetikJo commented Jul 27, 2023

FYI
A new api execution test with the optional attribute fac5467 7d7a5fa was successful.

@bnavigator bnavigator merged commit 25fa7b0 into python-control:master Aug 24, 2023
51 checks passed
@bnavigator bnavigator added this to the 0.6.0 milestone Aug 26, 2023
@bnavigator bnavigator mentioned this pull request Sep 22, 2023
4 tasks
@KybernetikJo KybernetikJo deleted the bugfix_ab13bd branch January 6, 2024 20:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants