lapack.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717
  1. """
  2. Low-level LAPACK functions (:mod:`scipy.linalg.lapack`)
  3. =======================================================
  4. This module contains low-level functions from the LAPACK library.
  5. The `*gegv` family of routines have been removed from LAPACK 3.6.0
  6. and have been deprecated in SciPy 0.17.0. They will be removed in
  7. a future release.
  8. .. versionadded:: 0.12.0
  9. .. note::
  10. The common ``overwrite_<>`` option in many routines, allows the
  11. input arrays to be overwritten to avoid extra memory allocation.
  12. However this requires the array to satisfy two conditions
  13. which are memory order and the data type to match exactly the
  14. order and the type expected by the routine.
  15. As an example, if you pass a double precision float array to any
  16. ``S....`` routine which expects single precision arguments, f2py
  17. will create an intermediate array to match the argument types and
  18. overwriting will be performed on that intermediate array.
  19. Similarly, if a C-contiguous array is passed, f2py will pass a
  20. FORTRAN-contiguous array internally. Please make sure that these
  21. details are satisfied. More information can be found in the f2py
  22. documentation.
  23. .. warning::
  24. These functions do little to no error checking.
  25. It is possible to cause crashes by mis-using them,
  26. so prefer using the higher-level routines in `scipy.linalg`.
  27. Finding functions
  28. -----------------
  29. .. autosummary::
  30. get_lapack_funcs
  31. All functions
  32. -------------
  33. .. autosummary::
  34. :toctree: generated/
  35. sgbsv
  36. dgbsv
  37. cgbsv
  38. zgbsv
  39. sgbtrf
  40. dgbtrf
  41. cgbtrf
  42. zgbtrf
  43. sgbtrs
  44. dgbtrs
  45. cgbtrs
  46. zgbtrs
  47. sgebal
  48. dgebal
  49. cgebal
  50. zgebal
  51. sgees
  52. dgees
  53. cgees
  54. zgees
  55. sgeev
  56. dgeev
  57. cgeev
  58. zgeev
  59. sgeev_lwork
  60. dgeev_lwork
  61. cgeev_lwork
  62. zgeev_lwork
  63. sgegv
  64. dgegv
  65. cgegv
  66. zgegv
  67. sgehrd
  68. dgehrd
  69. cgehrd
  70. zgehrd
  71. sgehrd_lwork
  72. dgehrd_lwork
  73. cgehrd_lwork
  74. zgehrd_lwork
  75. sgelss
  76. dgelss
  77. cgelss
  78. zgelss
  79. sgelss_lwork
  80. dgelss_lwork
  81. cgelss_lwork
  82. zgelss_lwork
  83. sgelsd
  84. dgelsd
  85. cgelsd
  86. zgelsd
  87. sgelsd_lwork
  88. dgelsd_lwork
  89. cgelsd_lwork
  90. zgelsd_lwork
  91. sgelsy
  92. dgelsy
  93. cgelsy
  94. zgelsy
  95. sgelsy_lwork
  96. dgelsy_lwork
  97. cgelsy_lwork
  98. zgelsy_lwork
  99. sgeqp3
  100. dgeqp3
  101. cgeqp3
  102. zgeqp3
  103. sgeqrf
  104. dgeqrf
  105. cgeqrf
  106. zgeqrf
  107. sgerqf
  108. dgerqf
  109. cgerqf
  110. zgerqf
  111. sgesdd
  112. dgesdd
  113. cgesdd
  114. zgesdd
  115. sgesdd_lwork
  116. dgesdd_lwork
  117. cgesdd_lwork
  118. zgesdd_lwork
  119. sgesvd
  120. dgesvd
  121. cgesvd
  122. zgesvd
  123. sgesvd_lwork
  124. dgesvd_lwork
  125. cgesvd_lwork
  126. zgesvd_lwork
  127. sgesv
  128. dgesv
  129. cgesv
  130. zgesv
  131. sgesvx
  132. dgesvx
  133. cgesvx
  134. zgesvx
  135. sgecon
  136. dgecon
  137. cgecon
  138. zgecon
  139. ssysv
  140. dsysv
  141. csysv
  142. zsysv
  143. ssysv_lwork
  144. dsysv_lwork
  145. csysv_lwork
  146. zsysv_lwork
  147. ssysvx
  148. dsysvx
  149. csysvx
  150. zsysvx
  151. ssysvx_lwork
  152. dsysvx_lwork
  153. csysvx_lwork
  154. zsysvx_lwork
  155. ssygst
  156. dsygst
  157. ssytrd
  158. dsytrd
  159. ssytrd_lwork
  160. dsytrd_lwork
  161. chetrd
  162. zhetrd
  163. chetrd_lwork
  164. zhetrd_lwork
  165. chesv
  166. zhesv
  167. chesv_lwork
  168. zhesv_lwork
  169. chesvx
  170. zhesvx
  171. chesvx_lwork
  172. zhesvx_lwork
  173. chegst
  174. zhegst
  175. sgetrf
  176. dgetrf
  177. cgetrf
  178. zgetrf
  179. sgetri
  180. dgetri
  181. cgetri
  182. zgetri
  183. sgetri_lwork
  184. dgetri_lwork
  185. cgetri_lwork
  186. zgetri_lwork
  187. sgetrs
  188. dgetrs
  189. cgetrs
  190. zgetrs
  191. sgges
  192. dgges
  193. cgges
  194. zgges
  195. sggev
  196. dggev
  197. cggev
  198. zggev
  199. chbevd
  200. zhbevd
  201. chbevx
  202. zhbevx
  203. cheev
  204. zheev
  205. cheevd
  206. zheevd
  207. cheevr
  208. zheevr
  209. chegv
  210. zhegv
  211. chegvd
  212. zhegvd
  213. chegvx
  214. zhegvx
  215. slarf
  216. dlarf
  217. clarf
  218. zlarf
  219. slarfg
  220. dlarfg
  221. clarfg
  222. zlarfg
  223. slartg
  224. dlartg
  225. clartg
  226. zlartg
  227. slasd4
  228. dlasd4
  229. slaswp
  230. dlaswp
  231. claswp
  232. zlaswp
  233. slauum
  234. dlauum
  235. clauum
  236. zlauum
  237. spbsv
  238. dpbsv
  239. cpbsv
  240. zpbsv
  241. spbtrf
  242. dpbtrf
  243. cpbtrf
  244. zpbtrf
  245. spbtrs
  246. dpbtrs
  247. cpbtrs
  248. zpbtrs
  249. sposv
  250. dposv
  251. cposv
  252. zposv
  253. sposvx
  254. dposvx
  255. cposvx
  256. zposvx
  257. spocon
  258. dpocon
  259. cpocon
  260. zpocon
  261. spotrf
  262. dpotrf
  263. cpotrf
  264. zpotrf
  265. spotri
  266. dpotri
  267. cpotri
  268. zpotri
  269. spotrs
  270. dpotrs
  271. cpotrs
  272. zpotrs
  273. crot
  274. zrot
  275. strsyl
  276. dtrsyl
  277. ctrsyl
  278. ztrsyl
  279. strtri
  280. dtrtri
  281. ctrtri
  282. ztrtri
  283. strtrs
  284. dtrtrs
  285. ctrtrs
  286. ztrtrs
  287. spftrf
  288. dpftrf
  289. cpftrf
  290. zpftrf
  291. spftri
  292. dpftri
  293. cpftri
  294. zpftri
  295. spftrs
  296. dpftrs
  297. cpftrs
  298. zpftrs
  299. cunghr
  300. zunghr
  301. cungqr
  302. zungqr
  303. cungrq
  304. zungrq
  305. cunmqr
  306. zunmqr
  307. cunmrz
  308. zunmrz
  309. cunmrz_lwork
  310. zunmrz_lwork
  311. sgtsv
  312. dgtsv
  313. cgtsv
  314. zgtsv
  315. sptsv
  316. dptsv
  317. cptsv
  318. zptsv
  319. slamch
  320. dlamch
  321. sorghr
  322. dorghr
  323. sorgqr
  324. dorgqr
  325. sorgrq
  326. dorgrq
  327. sormqr
  328. dormqr
  329. sormrz
  330. dormrz
  331. sormrz_lwork
  332. dormrz_lwork
  333. ssbev
  334. dsbev
  335. ssbevd
  336. dsbevd
  337. ssbevx
  338. dsbevx
  339. sstebz
  340. dstebz
  341. sstemr
  342. dstemr
  343. ssterf
  344. dsterf
  345. sstein
  346. dstein
  347. sstev
  348. dstev
  349. ssyev
  350. dsyev
  351. ssyevd
  352. dsyevd
  353. ssyevr
  354. dsyevr
  355. ssygv
  356. dsygv
  357. ssygvd
  358. dsygvd
  359. ssygvx
  360. dsygvx
  361. ssfrk
  362. dsfrk
  363. chfrk
  364. zhfrk
  365. stfsm
  366. dtfsm
  367. ctfsm
  368. ztfsm
  369. stpttf
  370. dtpttf
  371. ctpttf
  372. ztpttf
  373. stfttp
  374. dtfttp
  375. ctfttp
  376. ztfttp
  377. stfttr
  378. dtfttr
  379. ctfttr
  380. ztfttr
  381. strttf
  382. dtrttf
  383. ctrttf
  384. ztrttf
  385. stpttr
  386. dtpttr
  387. ctpttr
  388. ztpttr
  389. strttp
  390. dtrttp
  391. ctrttp
  392. ztrttp
  393. stfsm
  394. dtfsm
  395. ctfsm
  396. dtfsm
  397. stzrzf
  398. dtzrzf
  399. ctzrzf
  400. ztzrzf
  401. stzrzf_lwork
  402. dtzrzf_lwork
  403. ctzrzf_lwork
  404. ztzrzf_lwork
  405. slange
  406. dlange
  407. clange
  408. zlange
  409. ilaver
  410. """
  411. #
  412. # Author: Pearu Peterson, March 2002
  413. #
  414. from __future__ import division, print_function, absolute_import
  415. __all__ = ['get_lapack_funcs']
  416. import numpy as _np
  417. from .blas import _get_funcs
  418. # Backward compatibility:
  419. from .blas import find_best_blas_type as find_best_lapack_type
  420. from scipy.linalg import _flapack
  421. try:
  422. from scipy.linalg import _clapack
  423. except ImportError:
  424. _clapack = None
  425. # Backward compatibility
  426. from scipy._lib._util import DeprecatedImport as _DeprecatedImport
  427. clapack = _DeprecatedImport("scipy.linalg.blas.clapack", "scipy.linalg.lapack")
  428. flapack = _DeprecatedImport("scipy.linalg.blas.flapack", "scipy.linalg.lapack")
  429. # Expose all functions (only flapack --- clapack is an implementation detail)
  430. empty_module = None
  431. from scipy.linalg._flapack import *
  432. del empty_module
  433. _dep_message = """The `*gegv` family of routines has been deprecated in
  434. LAPACK 3.6.0 in favor of the `*ggev` family of routines.
  435. The corresponding wrappers will be removed from SciPy in
  436. a future release."""
  437. cgegv = _np.deprecate(cgegv, old_name='cgegv', message=_dep_message)
  438. dgegv = _np.deprecate(dgegv, old_name='dgegv', message=_dep_message)
  439. sgegv = _np.deprecate(sgegv, old_name='sgegv', message=_dep_message)
  440. zgegv = _np.deprecate(zgegv, old_name='zgegv', message=_dep_message)
  441. # Modyfy _flapack in this scope so the deprecation warnings apply to
  442. # functions returned by get_lapack_funcs.
  443. _flapack.cgegv = cgegv
  444. _flapack.dgegv = dgegv
  445. _flapack.sgegv = sgegv
  446. _flapack.zgegv = zgegv
  447. # some convenience alias for complex functions
  448. _lapack_alias = {
  449. 'corghr': 'cunghr', 'zorghr': 'zunghr',
  450. 'corghr_lwork': 'cunghr_lwork', 'zorghr_lwork': 'zunghr_lwork',
  451. 'corgqr': 'cungqr', 'zorgqr': 'zungqr',
  452. 'cormqr': 'cunmqr', 'zormqr': 'zunmqr',
  453. 'corgrq': 'cungrq', 'zorgrq': 'zungrq',
  454. }
  455. def get_lapack_funcs(names, arrays=(), dtype=None):
  456. """Return available LAPACK function objects from names.
  457. Arrays are used to determine the optimal prefix of LAPACK routines.
  458. Parameters
  459. ----------
  460. names : str or sequence of str
  461. Name(s) of LAPACK functions without type prefix.
  462. arrays : sequence of ndarrays, optional
  463. Arrays can be given to determine optimal prefix of LAPACK
  464. routines. If not given, double-precision routines will be
  465. used, otherwise the most generic type in arrays will be used.
  466. dtype : str or dtype, optional
  467. Data-type specifier. Not used if `arrays` is non-empty.
  468. Returns
  469. -------
  470. funcs : list
  471. List containing the found function(s).
  472. Notes
  473. -----
  474. This routine automatically chooses between Fortran/C
  475. interfaces. Fortran code is used whenever possible for arrays with
  476. column major order. In all other cases, C code is preferred.
  477. In LAPACK, the naming convention is that all functions start with a
  478. type prefix, which depends on the type of the principal
  479. matrix. These can be one of {'s', 'd', 'c', 'z'} for the numpy
  480. types {float32, float64, complex64, complex128} respectively, and
  481. are stored in attribute ``typecode`` of the returned functions.
  482. Examples
  483. --------
  484. Suppose we would like to use '?lange' routine which computes the selected
  485. norm of an array. We pass our array in order to get the correct 'lange'
  486. flavor.
  487. >>> import scipy.linalg as LA
  488. >>> a = np.random.rand(3,2)
  489. >>> x_lange = LA.get_lapack_funcs('lange', (a,))
  490. >>> x_lange.typecode
  491. 'd'
  492. >>> x_lange = LA.get_lapack_funcs('lange',(a*1j,))
  493. >>> x_lange.typecode
  494. 'z'
  495. Several LAPACK routines work best when its internal WORK array has
  496. the optimal size (big enough for fast computation and small enough to
  497. avoid waste of memory). This size is determined also by a dedicated query
  498. to the function which is often wrapped as a standalone function and
  499. commonly denoted as ``###_lwork``. Below is an example for ``?sysv``
  500. >>> import scipy.linalg as LA
  501. >>> a = np.random.rand(1000,1000)
  502. >>> b = np.random.rand(1000,1)*1j
  503. >>> # We pick up zsysv and zsysv_lwork due to b array
  504. ... xsysv, xlwork = LA.get_lapack_funcs(('sysv', 'sysv_lwork'), (a, b))
  505. >>> opt_lwork, _ = xlwork(a.shape[0]) # returns a complex for 'z' prefix
  506. >>> udut, ipiv, x, info = xsysv(a, b, lwork=int(opt_lwork.real))
  507. """
  508. return _get_funcs(names, arrays, dtype,
  509. "LAPACK", _flapack, _clapack,
  510. "flapack", "clapack", _lapack_alias)
  511. def _compute_lwork(routine, *args, **kwargs):
  512. """
  513. Round floating-point lwork returned by lapack to integer.
  514. Several LAPACK routines compute optimal values for LWORK, which
  515. they return in a floating-point variable. However, for large
  516. values of LWORK, single-precision floating point is not sufficient
  517. to hold the exact value --- some LAPACK versions (<= 3.5.0 at
  518. least) truncate the returned integer to single precision and in
  519. some cases this can be smaller than the required value.
  520. Examples
  521. --------
  522. >>> from scipy.linalg import lapack
  523. >>> n = 5000
  524. >>> s_r, s_lw = lapack.get_lapack_funcs(('sysvx', 'sysvx_lwork'))
  525. >>> lwork = lapack._compute_lwork(s_lw, n)
  526. >>> lwork
  527. 32000
  528. """
  529. wi = routine(*args, **kwargs)
  530. if len(wi) < 2:
  531. raise ValueError('')
  532. info = wi[-1]
  533. if info != 0:
  534. raise ValueError("Internal work array size computation failed: "
  535. "%d" % (info,))
  536. lwork = [w.real for w in wi[:-1]]
  537. dtype = getattr(routine, 'dtype', None)
  538. if dtype == _np.float32 or dtype == _np.complex64:
  539. # Single-precision routine -- take next fp value to work
  540. # around possible truncation in LAPACK code
  541. lwork = _np.nextafter(lwork, _np.inf, dtype=_np.float32)
  542. lwork = _np.array(lwork, _np.int64)
  543. if _np.any(_np.logical_or(lwork < 0, lwork > _np.iinfo(_np.int32).max)):
  544. raise ValueError("Too large work array required -- computation cannot "
  545. "be performed with standard 32-bit LAPACK.")
  546. lwork = lwork.astype(_np.int32)
  547. if lwork.size == 1:
  548. return lwork[0]
  549. return lwork