You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

blas_lapack.pyf.src 16 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417
  1. !
  2. ! Taken from scipy/linalg
  3. !
  4. ! Shorthand notations
  5. !
  6. ! <tchar=s,d,cs,zd>
  7. ! <tchar2c=cs,zd>
  8. !
  9. ! <prefix2=s,d>
  10. ! <prefix2c=c,z>
  11. ! <prefix3=s,sc>
  12. ! <prefix4=d,dz>
  13. ! <prefix6=s,d,c,z,c,z>
  14. !
  15. ! <ftype2=real,double precision>
  16. ! <ftype2c=complex,double complex>
  17. ! <ftype3=real,complex>
  18. ! <ftype4=double precision,double complex>
  19. ! <ftypereal3=real,real>
  20. ! <ftypereal4=double precision,double precision>
  21. ! <ftype6=real,double precision,complex,double complex,\2,\3>
  22. ! <ftype6creal=real,double precision,complex,double complex,\0,\1>
  23. !
  24. ! <ctype2=float,double>
  25. ! <ctype2c=complex_float,complex_double>
  26. ! <ctype3=float,complex_float>
  27. ! <ctype4=double,complex_double>
  28. ! <ctypereal3=float,float>
  29. ! <ctypereal4=double,double>
  30. ! <ctype6=float,double,complex_float,complex_double,\2,\3>
  31. ! <ctype6creal=float,double,complex_float,complex_double,\0,\1>
  32. !
  33. !
  34. ! Level 1 BLAS
  35. !
  36. python module _flapack
  37. usercode '''
  38. #define F_INT int
  39. '''
  40. interface
  41. subroutine <prefix>axpy(n,a,x,offx,incx,y,offy,incy)
  42. ! Calculate z = a*x+y, where a is scalar.
  43. callstatement (*f2py_func)(&n,&a,x+offx,&incx,y+offy,&incy)
  44. callprotoargument F_INT*,<ctype>*,<ctype>*,F_INT*,<ctype>*,F_INT*
  45. <ftype> dimension(*), intent(in) :: x
  46. <ftype> dimension(*), intent(in,out,out=z) :: y
  47. <ftype> optional, intent(in):: a=<1.0,\0,(1.0\,0.0),\2>
  48. integer optional, intent(in),check(incx>0||incx<0) :: incx = 1
  49. integer optional, intent(in),check(incy>0||incy<0) :: incy = 1
  50. integer optional, intent(in),depend(x) :: offx=0
  51. integer optional, intent(in),depend(y) :: offy=0
  52. check(offx>=0 && offx<len(x)) :: offx
  53. check(offy>=0 && offy<len(y)) :: offy
  54. integer optional, intent(in),depend(x,incx,offx,y,incy,offy) :: &
  55. n = (len(x)-offx)/abs(incx)
  56. check(len(x)-offx>(n-1)*abs(incx)) :: n
  57. check(len(y)-offy>(n-1)*abs(incy)) :: n
  58. end subroutine <prefix>axpy
  59. function ddot(n,x,offx,incx,y,offy,incy) result (xy)
  60. ! Computes a vector-vector dot product.
  61. callstatement ddot_return_value = (*f2py_func)(&n,x+offx,&incx,y+offy,&incy)
  62. callprotoargument F_INT*,double*,F_INT*,double*,F_INT*
  63. intent(c) ddot
  64. fortranname F_FUNC(ddot,DDOT)
  65. double precision dimension(*), intent(in) :: x
  66. double precision dimension(*), intent(in) :: y
  67. double precision ddot,xy
  68. integer optional, intent(in),check(incx>0||incx<0) :: incx = 1
  69. integer optional, intent(in),check(incy>0||incy<0) :: incy = 1
  70. integer optional, intent(in),depend(x) :: offx=0
  71. integer optional, intent(in),depend(y) :: offy=0
  72. check(offx>=0 && offx<len(x)) :: offx
  73. check(offy>=0 && offy<len(y)) :: offy
  74. integer optional, intent(in),depend(x,incx,offx,y,incy,offy) :: &
  75. n = (len(x)-offx)/abs(incx)
  76. check(len(x)-offx>(n-1)*abs(incx)) :: n
  77. check(len(y)-offy>(n-1)*abs(incy)) :: n
  78. end function ddot
  79. function <prefix4>nrm2(n,x,offx,incx) result(n2)
  80. <ftypereal4> <prefix4>nrm2, n2
  81. callstatement <prefix4>nrm2_return_value = (*f2py_func)(&n,x+offx,&incx)
  82. callprotoargument F_INT*,<ctype4>*,F_INT*
  83. intent(c) <prefix4>nrm2
  84. fortranname F_FUNC(<prefix4>nrm2,<D,DZ>NRM2)
  85. <ftype4> dimension(*),intent(in) :: x
  86. integer optional, intent(in),check(incx>0) :: incx = 1
  87. integer optional,intent(in),depend(x) :: offx=0
  88. check(offx>=0 && offx<len(x)) :: offx
  89. integer optional,intent(in),depend(x,incx,offx) :: n = (len(x)-offx)/abs(incx)
  90. check(len(x)-offx>(n-1)*abs(incx)) :: n
  91. end function <prefix4>nrm2
  92. !
  93. ! Level 2 BLAS
  94. !
  95. subroutine <prefix>gemv(m,n,alpha,a,x,beta,y,offx,incx,offy,incy,trans,rows,cols,ly)
  96. ! Computes a matrix-vector product using a general matrix
  97. !
  98. ! y = gemv(alpha,a,x,beta=0,y=0,offx=0,incx=1,offy=0,incy=0,trans=0)
  99. ! Calculate y <- alpha * op(A) * x + beta * y
  100. callstatement (*f2py_func)((trans?(trans==2?"C":"T"):"N"),&m,&n,&alpha,a,&m, &
  101. x+offx,&incx,&beta,y+offy,&incy)
  102. callprotoargument char*,F_INT*,F_INT*,<ctype>*,<ctype>*,F_INT*,<ctype>*,F_INT*,<ctype>*, &
  103. <ctype>*,F_INT*
  104. integer optional, intent(in), check(trans>=0 && trans <=2) :: trans = 0
  105. integer optional, intent(in), check(incx>0||incx<0) :: incx = 1
  106. integer optional, intent(in), check(incy>0||incy<0) :: incy = 1
  107. <ftype> intent(in) :: alpha
  108. <ftype> intent(in), optional :: beta = <0.0,\0,(0.0\,0.0),\2>
  109. <ftype> dimension(*), intent(in) :: x
  110. <ftype> dimension(ly), intent(in,copy,out), depend(ly),optional :: y
  111. integer intent(hide), depend(incy,rows,offy) :: ly = &
  112. (y_capi==Py_None?1+offy+(rows-1)*abs(incy):-1)
  113. <ftype> dimension(m,n), intent(in) :: a
  114. integer depend(a), intent(hide):: m = shape(a,0)
  115. integer depend(a), intent(hide):: n = shape(a,1)
  116. integer optional, intent(in) :: offx=0
  117. integer optional, intent(in) :: offy=0
  118. check(offx>=0 && offx<len(x)) :: x
  119. check(len(x)>offx+(cols-1)*abs(incx)) :: x
  120. depend(offx,cols,incx) :: x
  121. check(offy>=0 && offy<len(y)) :: y
  122. check(len(y)>offy+(rows-1)*abs(incy)) :: y
  123. depend(offy,rows,incy) :: y
  124. integer depend(m,n,trans), intent(hide) :: rows = (trans?n:m)
  125. integer depend(m,n,trans), intent(hide) :: cols = (trans?m:n)
  126. end subroutine <prefix>gemv
  127. subroutine <prefix>gbmv(m,n,kl,ku,alpha,a,lda,x,incx,offx,beta,y,incy,offy,trans,ly)
  128. ! Performs one of the matrix-vector operations
  129. !
  130. ! y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y,
  131. ! or y := alpha*A**H*x + beta*y,
  132. !
  133. ! where alpha and beta are scalars, x and y are vectors and A is an
  134. ! m by n band matrix, with kl sub-diagonals and ku super-diagonals.
  135. callstatement (*f2py_func)((trans?(trans==2?"C":"T"):"N"),&m,&n,&kl,&ku,&alpha,a,&lda,x+offx,&incx,&beta,y+offy,&incy)
  136. callprotoargument char*,F_INT*,F_INT*,F_INT*,F_INT*,<ctype>*,<ctype>*,F_INT*,<ctype>*,F_INT*,<ctype>*,<ctype>*,F_INT*
  137. integer optional,intent(in),check(trans>=0 && trans <=2) :: trans = 0
  138. integer intent(in), depend(ku,kl),check(m>=ku+kl+1) :: m
  139. integer intent(in),check(n>=0&&n==shape(a,1)),depend(a) :: n
  140. integer intent(in),check(kl>=0) :: kl
  141. integer intent(in),check(ku>=0) :: ku
  142. integer intent(hide),depend(a) :: lda = MAX(shape(a,0),1)
  143. integer optional, intent(in),check(incx>0||incx<0) :: incx = 1
  144. integer optional, intent(in),check(incy>0||incy<0) :: incy = 1
  145. integer intent(hide),depend(m,n,incy,offy,trans) :: ly = &
  146. (y_capi==Py_None?1+offy+(trans==0?m-1:n-1)*abs(incy):-1)
  147. integer optional, intent(in) :: offx=0
  148. integer optional, intent(in) :: offy=0
  149. <ftype> intent(in) :: alpha
  150. <ftype> intent(in),optional :: beta = <0.0,\0,(0.0\,0.0),\2>
  151. <ftype> dimension(lda,n),intent(in) :: a
  152. <ftype> dimension(ly), intent(in,out,copy,out=yout),depend(ly),optional :: y
  153. check(offy>=0 && offy<len(y)) :: y
  154. check(len(y)>offy+(trans==0?m-1:n-1)*abs(incy)) :: y
  155. depend(offy,n,incy) :: y
  156. <ftype> dimension(*), intent(in) :: x
  157. check(offx>=0 && offx<len(x)) :: x
  158. check(len(x)>offx+(trans==0?n-1:m-1)*abs(incx)) :: x
  159. depend(offx,n,incx) :: x
  160. end subroutine <prefix>gbmv
  161. !
  162. ! Level 3 BLAS
  163. !
  164. subroutine <prefix>gemm(m,n,k,alpha,a,b,beta,c,trans_a,trans_b,lda,ka,ldb,kb)
  165. ! Computes a scalar-matrix-matrix product and adds the result to a
  166. ! scalar-matrix product.
  167. !
  168. ! c = gemm(alpha,a,b,beta=0,c=0,trans_a=0,trans_b=0,overwrite_c=0)
  169. ! Calculate C <- alpha * op(A) * op(B) + beta * C
  170. callstatement (*f2py_func)((trans_a?(trans_a==2?"C":"T"):"N"), &
  171. (trans_b?(trans_b==2?"C":"T"):"N"),&m,&n,&k,&alpha,a,&lda,b,&ldb,&beta,c,&m)
  172. callprotoargument char*,char*,F_INT*,F_INT*,F_INT*,<ctype>*,<ctype>*,F_INT*,<ctype>*, &
  173. F_INT*,<ctype>*,<ctype>*,F_INT*
  174. integer optional,intent(in),check(trans_a>=0 && trans_a <=2) :: trans_a = 0
  175. integer optional,intent(in),check(trans_b>=0 && trans_b <=2) :: trans_b = 0
  176. <ftype> intent(in) :: alpha
  177. <ftype> intent(in),optional :: beta = <0.0,\0,(0.0\,0.0),\2>
  178. <ftype> dimension(lda,ka),intent(in) :: a
  179. <ftype> dimension(ldb,kb),intent(in) :: b
  180. <ftype> dimension(m,n),intent(in,out,copy),depend(m,n),optional :: c
  181. check(shape(c,0)==m && shape(c,1)==n) :: c
  182. integer depend(a),intent(hide) :: lda = shape(a,0)
  183. integer depend(a),intent(hide) :: ka = shape(a,1)
  184. integer depend(b),intent(hide) :: ldb = shape(b,0)
  185. integer depend(b),intent(hide) :: kb = shape(b,1)
  186. integer depend(a,trans_a,ka,lda),intent(hide):: m = (trans_a?ka:lda)
  187. integer depend(a,trans_a,ka,lda),intent(hide):: k = (trans_a?lda:ka)
  188. integer depend(b,trans_b,kb,ldb,k),intent(hide),check(trans_b?kb==k:ldb==k) :: &
  189. n = (trans_b?ldb:kb)
  190. end subroutine <prefix>gemm
  191. subroutine <prefix6><sy,\0,\0,\0,he,he>rk(n,k,alpha,a,beta,c,trans,lower,lda,ka)
  192. ! performs one of the symmetric rank k operations
  193. ! C := alpha*A*A**T + beta*C, or C := alpha*A**T*A + beta*C,
  194. !
  195. ! c = syrk(alpha,a,beta=0,c=0,trans=0,lower=0,overwrite_c=0)
  196. !
  197. callstatement (*f2py_func)((lower?"L":"U"), &
  198. (trans?(trans==2?"C":"T"):"N"), &n,&k,&alpha,a,&lda,&beta,c,&n)
  199. callprotoargument char*,char*,F_INT*,F_INT*,<ctype6>*,<ctype6>*,F_INT*,<ctype6>*, &
  200. <ctype6>*,F_INT*
  201. integer optional, intent(in),check(lower==0||lower==1) :: lower = 0
  202. integer optional,intent(in),check(trans>=0 && trans <=2) :: trans = 0
  203. <ftype6> intent(in) :: alpha
  204. <ftype6> intent(in),optional :: beta = <0.0,\0,(0.0\,0.0),\2,\2,\2>
  205. <ftype6> dimension(lda,ka),intent(in) :: a
  206. <ftype6> dimension(n,n),intent(in,out,copy),depend(n),optional :: c
  207. check(shape(c,0)==n && shape(c,1)==n) :: c
  208. integer depend(a),intent(hide) :: lda = shape(a,0)
  209. integer depend(a),intent(hide) :: ka = shape(a,1)
  210. integer depend(a, trans, ka, lda), intent(hide) :: n = (trans ? ka : lda)
  211. integer depend(a, trans, ka, lda), intent(hide) :: k = (trans ? lda : ka)
  212. end subroutine <prefix6><sy,\0,\0,\0,he,he>rk
  213. !
  214. ! LAPACK
  215. !
  216. subroutine <prefix>gesv(n,nrhs,a,piv,b,info)
  217. ! lu,piv,x,info = gesv(a,b,overwrite_a=0,overwrite_b=0)
  218. ! Solve A * X = B.
  219. ! A = P * L * U
  220. ! U is upper diagonal triangular, L is unit lower triangular,
  221. ! piv pivots columns.
  222. callstatement {F_INT i;(*f2py_func)(&n,&nrhs,a,&n,piv,b,&n,&info);for(i=0;i\<n;--piv[i++]);}
  223. callprotoargument F_INT*,F_INT*,<ctype>*,F_INT*,F_INT*,<ctype>*,F_INT*,F_INT*
  224. integer depend(a),intent(hide):: n = shape(a,0)
  225. integer depend(b),intent(hide):: nrhs = shape(b,1)
  226. <ftype> dimension(n,n),check(shape(a,0)==shape(a,1)) :: a
  227. integer dimension(n),depend(n),intent(out) :: piv
  228. <ftype> dimension(n,nrhs),check(shape(a,0)==shape(b,0)),depend(n) :: b
  229. integer intent(out)::info
  230. intent(in,out,copy,out=x) b
  231. intent(in,out,copy,out=lu) a
  232. end subroutine <prefix>gesv
  233. subroutine <prefix2>gesdd(m,n,minmn,u0,u1,vt0,vt1,a,compute_uv,full_matrices,u,s,vt,work,lwork,iwork,info)
  234. ! u,s,vt,info = gesdd(a,compute_uv=1,lwork=..,overwrite_a=0)
  235. ! Compute the singular value decomposition (SVD) using divide and conquer:
  236. ! A = U * SIGMA * transpose(V)
  237. ! A - M x N matrix
  238. ! U - M x M matrix or min(M,N) x N if full_matrices=False
  239. ! SIGMA - M x N zero matrix with a main diagonal filled with min(M,N)
  240. ! singular values
  241. ! transpose(V) - N x N matrix or N x min(M,N) if full_matrices=False
  242. callstatement (*f2py_func)((compute_uv?(full_matrices?"A":"S"):"N"),&m,&n,a,&m,s,u,&u0,vt,&vt0,work,&lwork,iwork,&info)
  243. callprotoargument char*,F_INT*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,<ctype2>*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,F_INT*,F_INT*,F_INT*
  244. integer intent(in),optional,check(compute_uv==0||compute_uv==1):: compute_uv = 1
  245. integer intent(in),optional,check(full_matrices==0||full_matrices==1):: full_matrices = 1
  246. integer intent(hide),depend(a):: m = shape(a,0)
  247. integer intent(hide),depend(a):: n = shape(a,1)
  248. integer intent(hide),depend(m,n):: minmn = MIN(m,n)
  249. integer intent(hide),depend(compute_uv,minmn) :: u0 = (compute_uv?m:1)
  250. integer intent(hide),depend(compute_uv,minmn, full_matrices) :: u1 = (compute_uv?(full_matrices?m:minmn):1)
  251. integer intent(hide),depend(compute_uv,minmn, full_matrices) :: vt0 = (compute_uv?(full_matrices?n:minmn):1)
  252. integer intent(hide),depend(compute_uv,minmn) :: vt1 = (compute_uv?n:1)
  253. <ftype2> dimension(m,n),intent(in,copy,aligned8) :: a
  254. <ftype2> dimension(minmn),intent(out),depend(minmn) :: s
  255. <ftype2> dimension(u0,u1),intent(out),depend(u0, u1) :: u
  256. <ftype2> dimension(vt0,vt1),intent(out),depend(vt0, vt1) :: vt
  257. <ftype2> dimension(lwork),intent(hide,cache),depend(lwork) :: work
  258. integer optional,intent(in),depend(minmn,compute_uv) &
  259. :: lwork = max((compute_uv?4*minmn*minmn+MAX(m,n)+9*minmn:MAX(14*minmn+4,10*minmn+2+25*(25+8))+MAX(m,n)),1)
  260. integer intent(hide,cache),dimension(8*minmn),depend(minmn) :: iwork
  261. integer intent(out)::info
  262. end subroutine <prefix2>gesdd
  263. subroutine <prefix2>gesdd_lwork(m,n,minmn,u0,vt0,a,compute_uv,full_matrices,u,s,vt,work,lwork,iwork,info)
  264. ! LWORK computation for (S/D)GESDD
  265. fortranname <prefix2>gesdd
  266. callstatement (*f2py_func)((compute_uv?(full_matrices?"A":"S"):"N"),&m,&n,&a,&m,&s,&u,&u0,&vt,&vt0,&work,&lwork,&iwork,&info)
  267. callprotoargument char*,F_INT*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,<ctype2>*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,F_INT*,F_INT*,F_INT*
  268. integer intent(in),optional,check(compute_uv==0||compute_uv==1):: compute_uv = 1
  269. integer intent(in),optional,check(full_matrices==0||full_matrices==1):: full_matrices = 1
  270. integer intent(in) :: m
  271. integer intent(in) :: n
  272. integer intent(hide),depend(m,n):: minmn = MIN(m,n)
  273. integer intent(hide),depend(compute_uv,minmn) :: u0 = (compute_uv?m:1)
  274. integer intent(hide),depend(compute_uv,minmn, full_matrices) :: vt0 = (compute_uv?(full_matrices?n:minmn):1)
  275. <ftype2> intent(hide) :: a
  276. <ftype2> intent(hide) :: s
  277. <ftype2> intent(hide) :: u
  278. <ftype2> intent(hide) :: vt
  279. <ftype2> intent(out) :: work
  280. integer intent(hide) :: lwork = -1
  281. integer intent(hide) :: iwork
  282. integer intent(out) :: info
  283. end subroutine <prefix2>gesdd_lwork
  284. subroutine <prefix2>syev(compute_v,lower,n,w,a,lda,work,lwork,info)
  285. ! w,v,info = syev(a,compute_v=1,lower=0,lwork=3*n-1,overwrite_a=0)
  286. ! Compute all eigenvalues and, optionally, eigenvectors of a
  287. ! real symmetric matrix A.
  288. !
  289. ! Performance tip:
  290. ! If compute_v=0 then set also overwrite_a=1.
  291. callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,a,&lda,w,work,&lwork,&info)
  292. callprotoargument char*,char*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,<ctype2>*,F_INT*,F_INT*
  293. integer optional,intent(in):: compute_v = 1
  294. check(compute_v==1||compute_v==0) compute_v
  295. integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
  296. integer intent(hide),depend(a):: n = shape(a,0)
  297. integer intent(hide),depend(a):: lda = MAX(1,shape(a,0))
  298. <ftype2> dimension(n,n),check(shape(a,0)==shape(a,1)) :: a
  299. intent(in,copy,out,out=v) :: a
  300. <ftype2> dimension(n),intent(out),depend(n) :: w
  301. integer optional,intent(in),depend(n) :: lwork=max(3*n-1,1)
  302. check(lwork>=3*n-1) :: lwork
  303. <ftype2> dimension(lwork),intent(hide),depend(lwork) :: work
  304. integer intent(out) :: info
  305. end subroutine <prefix2>syev
  306. subroutine <prefix2>syev_lwork(lower,n,w,a,lda,work,lwork,info)
  307. ! LWORK routines for syev
  308. fortranname <prefix2>syev
  309. callstatement (*f2py_func)("N",(lower?"L":"U"),&n,&a,&lda,&w,&work,&lwork,&info)
  310. callprotoargument char*,char*,F_INT*,<ctype2>*,F_INT*,<ctype2>*,<ctype2>*,F_INT*,F_INT*
  311. integer intent(in):: n
  312. integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
  313. integer intent(hide),depend(n):: lda = MAX(1, n)
  314. <ftype2> intent(hide):: a
  315. <ftype2> intent(hide):: w
  316. integer intent(hide):: lwork = -1
  317. <ftype2> intent(out):: work
  318. integer intent(out):: info
  319. end subroutine <prefix2>syev_lwork
  320. end interface
  321. end python module _flapack