-
-
Notifications
You must be signed in to change notification settings - Fork 5.1k
/
lapack.py
1036 lines (797 loc) · 15.3 KB
/
lapack.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
Low-level LAPACK functions (:mod:`scipy.linalg.lapack`)
=======================================================
This module contains low-level functions from the LAPACK library.
.. versionadded:: 0.12.0
.. note::
The common ``overwrite_<>`` option in many routines, allows the
input arrays to be overwritten to avoid extra memory allocation.
However this requires the array to satisfy two conditions
which are memory order and the data type to match exactly the
order and the type expected by the routine.
As an example, if you pass a double precision float array to any
``S....`` routine which expects single precision arguments, f2py
will create an intermediate array to match the argument types and
overwriting will be performed on that intermediate array.
Similarly, if a C-contiguous array is passed, f2py will pass a
FORTRAN-contiguous array internally. Please make sure that these
details are satisfied. More information can be found in the f2py
documentation.
.. warning::
These functions do little to no error checking.
It is possible to cause crashes by mis-using them,
so prefer using the higher-level routines in `scipy.linalg`.
Finding functions
-----------------
.. autosummary::
:toctree: generated/
get_lapack_funcs
All functions
-------------
.. autosummary::
:toctree: generated/
sgbsv
dgbsv
cgbsv
zgbsv
sgbtrf
dgbtrf
cgbtrf
zgbtrf
sgbtrs
dgbtrs
cgbtrs
zgbtrs
sgebal
dgebal
cgebal
zgebal
sgecon
dgecon
cgecon
zgecon
sgeequ
dgeequ
cgeequ
zgeequ
sgeequb
dgeequb
cgeequb
zgeequb
sgees
dgees
cgees
zgees
sgeev
dgeev
cgeev
zgeev
sgeev_lwork
dgeev_lwork
cgeev_lwork
zgeev_lwork
sgehrd
dgehrd
cgehrd
zgehrd
sgehrd_lwork
dgehrd_lwork
cgehrd_lwork
zgehrd_lwork
sgejsv
dgejsv
sgels
dgels
cgels
zgels
sgels_lwork
dgels_lwork
cgels_lwork
zgels_lwork
sgelsd
dgelsd
cgelsd
zgelsd
sgelsd_lwork
dgelsd_lwork
cgelsd_lwork
zgelsd_lwork
sgelss
dgelss
cgelss
zgelss
sgelss_lwork
dgelss_lwork
cgelss_lwork
zgelss_lwork
sgelsy
dgelsy
cgelsy
zgelsy
sgelsy_lwork
dgelsy_lwork
cgelsy_lwork
zgelsy_lwork
sgeqp3
dgeqp3
cgeqp3
zgeqp3
sgeqrf
dgeqrf
cgeqrf
zgeqrf
sgeqrf_lwork
dgeqrf_lwork
cgeqrf_lwork
zgeqrf_lwork
sgeqrfp
dgeqrfp
cgeqrfp
zgeqrfp
sgeqrfp_lwork
dgeqrfp_lwork
cgeqrfp_lwork
zgeqrfp_lwork
sgerqf
dgerqf
cgerqf
zgerqf
sgesdd
dgesdd
cgesdd
zgesdd
sgesdd_lwork
dgesdd_lwork
cgesdd_lwork
zgesdd_lwork
sgesv
dgesv
cgesv
zgesv
sgesvd
dgesvd
cgesvd
zgesvd
sgesvd_lwork
dgesvd_lwork
cgesvd_lwork
zgesvd_lwork
sgesvx
dgesvx
cgesvx
zgesvx
sgetrf
dgetrf
cgetrf
zgetrf
sgetc2
dgetc2
cgetc2
zgetc2
sgetri
dgetri
cgetri
zgetri
sgetri_lwork
dgetri_lwork
cgetri_lwork
zgetri_lwork
sgetrs
dgetrs
cgetrs
zgetrs
sgesc2
dgesc2
cgesc2
zgesc2
sgges
dgges
cgges
zgges
sggev
dggev
cggev
zggev
sgglse
dgglse
cgglse
zgglse
sgglse_lwork
dgglse_lwork
cgglse_lwork
zgglse_lwork
sgtsv
dgtsv
cgtsv
zgtsv
sgtsvx
dgtsvx
cgtsvx
zgtsvx
chbevd
zhbevd
chbevx
zhbevx
checon
zhecon
cheequb
zheequb
cheev
zheev
cheev_lwork
zheev_lwork
cheevd
zheevd
cheevd_lwork
zheevd_lwork
cheevr
zheevr
cheevr_lwork
zheevr_lwork
cheevx
zheevx
cheevx_lwork
zheevx_lwork
chegst
zhegst
chegv
zhegv
chegv_lwork
zhegv_lwork
chegvd
zhegvd
chegvx
zhegvx
chegvx_lwork
zhegvx_lwork
chesv
zhesv
chesv_lwork
zhesv_lwork
chesvx
zhesvx
chesvx_lwork
zhesvx_lwork
chetrd
zhetrd
chetrd_lwork
zhetrd_lwork
chetrf
zhetrf
chetrf_lwork
zhetrf_lwork
chfrk
zhfrk
slamch
dlamch
slange
dlange
clange
zlange
slarf
dlarf
clarf
zlarf
slarfg
dlarfg
clarfg
zlarfg
slartg
dlartg
clartg
zlartg
slasd4
dlasd4
slaswp
dlaswp
claswp
zlaswp
slauum
dlauum
clauum
zlauum
sorcsd
dorcsd
sorcsd_lwork
dorcsd_lwork
sorghr
dorghr
sorghr_lwork
dorghr_lwork
sorgqr
dorgqr
sorgrq
dorgrq
sormqr
dormqr
sormrz
dormrz
sormrz_lwork
dormrz_lwork
spbsv
dpbsv
cpbsv
zpbsv
spbtrf
dpbtrf
cpbtrf
zpbtrf
spbtrs
dpbtrs
cpbtrs
zpbtrs
spftrf
dpftrf
cpftrf
zpftrf
spftri
dpftri
cpftri
zpftri
spftrs
dpftrs
cpftrs
zpftrs
spocon
dpocon
cpocon
zpocon
spstrf
dpstrf
cpstrf
zpstrf
spstf2
dpstf2
cpstf2
zpstf2
sposv
dposv
cposv
zposv
sposvx
dposvx
cposvx
zposvx
spotrf
dpotrf
cpotrf
zpotrf
spotri
dpotri
cpotri
zpotri
spotrs
dpotrs
cpotrs
zpotrs
sppcon
dppcon
cppcon
zppcon
sppsv
dppsv
cppsv
zppsv
spptrf
dpptrf
cpptrf
zpptrf
spptri
dpptri
cpptri
zpptri
spptrs
dpptrs
cpptrs
zpptrs
sptsv
dptsv
cptsv
zptsv
sptsvx
dptsvx
cptsvx
zptsvx
spttrf
dpttrf
cpttrf
zpttrf
spttrs
dpttrs
cpttrs
zpttrs
spteqr
dpteqr
cpteqr
zpteqr
crot
zrot
ssbev
dsbev
ssbevd
dsbevd
ssbevx
dsbevx
ssfrk
dsfrk
sstebz
dstebz
sstein
dstein
sstemr
dstemr
sstemr_lwork
dstemr_lwork
ssterf
dsterf
sstev
dstev
ssycon
dsycon
csycon
zsycon
ssyconv
dsyconv
csyconv
zsyconv
ssyequb
dsyequb
csyequb
zsyequb
ssyev
dsyev
ssyev_lwork
dsyev_lwork
ssyevd
dsyevd
ssyevd_lwork
dsyevd_lwork
ssyevr
dsyevr
ssyevr_lwork
dsyevr_lwork
ssyevx
dsyevx
ssyevx_lwork
dsyevx_lwork
ssygst
dsygst
ssygv
dsygv
ssygv_lwork
dsygv_lwork
ssygvd
dsygvd
ssygvx
dsygvx
ssygvx_lwork
dsygvx_lwork
ssysv
dsysv
csysv
zsysv
ssysv_lwork
dsysv_lwork
csysv_lwork
zsysv_lwork
ssysvx
dsysvx
csysvx
zsysvx
ssysvx_lwork
dsysvx_lwork
csysvx_lwork
zsysvx_lwork
ssytf2
dsytf2
csytf2
zsytf2
ssytrd
dsytrd
ssytrd_lwork
dsytrd_lwork
ssytrf
dsytrf
csytrf
zsytrf
ssytrf_lwork
dsytrf_lwork
csytrf_lwork
zsytrf_lwork
stbtrs
dtbtrs
ctbtrs
ztbtrs
stfsm
dtfsm
ctfsm
ztfsm
stfttp
dtfttp
ctfttp
ztfttp
stfttr
dtfttr
ctfttr
ztfttr
stgexc
dtgexc
ctgexc
ztgexc
stgsen
dtgsen
ctgsen
ztgsen
stgsen_lwork
dtgsen_lwork
ctgsen_lwork
ztgsen_lwork
stpttf
dtpttf
ctpttf
ztpttf
stpttr
dtpttr
ctpttr
ztpttr
strexc
dtrexc
ctrexc
ztrexc
strsen
dtrsen
ctrsen
ztrsen
strsen_lwork
dtrsen_lwork
ctrsen_lwork
ztrsen_lwork
strsyl
dtrsyl
ctrsyl
ztrsyl
strtri
dtrtri
ctrtri
ztrtri
strtrs
dtrtrs
ctrtrs
ztrtrs
strttf
dtrttf
ctrttf
ztrttf
strttp
dtrttp
ctrttp
ztrttp
stzrzf
dtzrzf
ctzrzf
ztzrzf
stzrzf_lwork
dtzrzf_lwork
ctzrzf_lwork
ztzrzf_lwork
cunghr
zunghr
cunghr_lwork
zunghr_lwork
cungqr
zungqr
cungrq
zungrq
cunmqr
zunmqr
sgeqrt
dgeqrt
cgeqrt
zgeqrt
sgemqrt
dgemqrt
cgemqrt
zgemqrt
sgttrf
dgttrf
cgttrf
zgttrf
sgttrs
dgttrs
cgttrs
zgttrs
stpqrt
dtpqrt
ctpqrt
ztpqrt
stpmqrt
dtpmqrt
ctpmqrt
ztpmqrt
cuncsd
zuncsd
cuncsd_lwork
zuncsd_lwork
cunmrz
zunmrz
cunmrz_lwork
zunmrz_lwork
ilaver
"""
#
# Author: Pearu Peterson, March 2002
#
import numpy as _np
from .blas import _get_funcs, _memoize_get_funcs
from scipy.linalg import _flapack
from re import compile as regex_compile
try:
from scipy.linalg import _clapack
except ImportError:
_clapack = None
try:
from scipy.linalg import _flapack_64
HAS_ILP64 = True
except ImportError:
HAS_ILP64 = False
_flapack_64 = None
# Expose all functions (only flapack --- clapack is an implementation detail)
empty_module = None
from scipy.linalg._flapack import *
del empty_module
__all__ = ['get_lapack_funcs']
# some convenience alias for complex functions
_lapack_alias = {
'corghr': 'cunghr', 'zorghr': 'zunghr',
'corghr_lwork': 'cunghr_lwork', 'zorghr_lwork': 'zunghr_lwork',
'corgqr': 'cungqr', 'zorgqr': 'zungqr',
'cormqr': 'cunmqr', 'zormqr': 'zunmqr',
'corgrq': 'cungrq', 'zorgrq': 'zungrq',
}
# Place guards against docstring rendering issues with special characters
p1 = regex_compile(r'with bounds (?P<b>.*?)( and (?P<s>.*?) storage){0,1}\n')
p2 = regex_compile(r'Default: (?P<d>.*?)\n')
def backtickrepl(m):
if m.group('s'):
return ('with bounds ``{}`` with ``{}`` storage\n'
''.format(m.group('b'), m.group('s')))
else:
return 'with bounds ``{}``\n'.format(m.group('b'))
for routine in [ssyevr, dsyevr, cheevr, zheevr,
ssyevx, dsyevx, cheevx, zheevx,
ssygvd, dsygvd, chegvd, zhegvd]:
if routine.__doc__:
routine.__doc__ = p1.sub(backtickrepl, routine.__doc__)
routine.__doc__ = p2.sub('Default ``\\1``\n', routine.__doc__)
else:
continue
del regex_compile, p1, p2, backtickrepl
@_memoize_get_funcs
def get_lapack_funcs(names, arrays=(), dtype=None, ilp64=False):
"""Return available LAPACK function objects from names.
Arrays are used to determine the optimal prefix of LAPACK routines.
Parameters
----------
names : str or sequence of str
Name(s) of LAPACK functions without type prefix.
arrays : sequence of ndarrays, optional
Arrays can be given to determine optimal prefix of LAPACK
routines. If not given, double-precision routines will be
used, otherwise the most generic type in arrays will be used.
dtype : str or dtype, optional
Data-type specifier. Not used if `arrays` is non-empty.
ilp64 : {True, False, 'preferred'}, optional
Whether to return ILP64 routine variant.
Choosing 'preferred' returns ILP64 routine if available, and
otherwise the 32-bit routine. Default: False
Returns
-------
funcs : list
List containing the found function(s).
Notes
-----
This routine automatically chooses between Fortran/C
interfaces. Fortran code is used whenever possible for arrays with
column major order. In all other cases, C code is preferred.
In LAPACK, the naming convention is that all functions start with a
type prefix, which depends on the type of the principal
matrix. These can be one of {'s', 'd', 'c', 'z'} for the NumPy
types {float32, float64, complex64, complex128} respectively, and
are stored in attribute ``typecode`` of the returned functions.
Examples
--------
Suppose we would like to use '?lange' routine which computes the selected
norm of an array. We pass our array in order to get the correct 'lange'
flavor.
>>> import numpy as np
>>> import scipy.linalg as LA
>>> rng = np.random.default_rng()
>>> a = rng.random((3,2))
>>> x_lange = LA.get_lapack_funcs('lange', (a,))
>>> x_lange.typecode
'd'
>>> x_lange = LA.get_lapack_funcs('lange',(a*1j,))
>>> x_lange.typecode
'z'
Several LAPACK routines work best when its internal WORK array has
the optimal size (big enough for fast computation and small enough to
avoid waste of memory). This size is determined also by a dedicated query
to the function which is often wrapped as a standalone function and
commonly denoted as ``###_lwork``. Below is an example for ``?sysv``
>>> a = rng.random((1000, 1000))
>>> b = rng.random((1000, 1)) * 1j
>>> # We pick up zsysv and zsysv_lwork due to b array
... xsysv, xlwork = LA.get_lapack_funcs(('sysv', 'sysv_lwork'), (a, b))
>>> opt_lwork, _ = xlwork(a.shape[0]) # returns a complex for 'z' prefix
>>> udut, ipiv, x, info = xsysv(a, b, lwork=int(opt_lwork.real))
"""
if isinstance(ilp64, str):
if ilp64 == 'preferred':
ilp64 = HAS_ILP64
else:
raise ValueError("Invalid value for 'ilp64'")
if not ilp64:
return _get_funcs(names, arrays, dtype,
"LAPACK", _flapack, _clapack,
"flapack", "clapack", _lapack_alias,
ilp64=False)
else:
if not HAS_ILP64:
raise RuntimeError("LAPACK ILP64 routine requested, but Scipy "
"compiled only with 32-bit BLAS")
return _get_funcs(names, arrays, dtype,
"LAPACK", _flapack_64, None,
"flapack_64", None, _lapack_alias,
ilp64=True)
_int32_max = _np.iinfo(_np.int32).max
_int64_max = _np.iinfo(_np.int64).max
def _compute_lwork(routine, *args, **kwargs):
"""
Round floating-point lwork returned by lapack to integer.
Several LAPACK routines compute optimal values for LWORK, which
they return in a floating-point variable. However, for large
values of LWORK, single-precision floating point is not sufficient
to hold the exact value --- some LAPACK versions (<= 3.5.0 at
least) truncate the returned integer to single precision and in
some cases this can be smaller than the required value.
Examples
--------
>>> from scipy.linalg import lapack
>>> n = 5000
>>> s_r, s_lw = lapack.get_lapack_funcs(('sysvx', 'sysvx_lwork'))
>>> lwork = lapack._compute_lwork(s_lw, n)
>>> lwork
32000
"""
dtype = getattr(routine, 'dtype', None)