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
|
## Automatically adapted for scipy Oct 18, 2005 by
#
# Author: Pearu Peterson, March 2002
#
__all__ = ['get_lapack_funcs']
import new
# The following ensures that possibly missing flavor (C or Fortran) is
# replaced with the available one. If none is available, exception
# is raised at the first attempt to use the resources.
import numpy
from scipy.linalg import flapack
from scipy.linalg import clapack
_use_force_clapack = 1
if hasattr(clapack,'empty_module'):
clapack = flapack
_use_force_clapack = 0
elif hasattr(flapack,'empty_module'):
flapack = clapack
def cast_to_lapack_prefix(t):
if issubclass(t, numpy.single):
prefix = 's'
elif issubclass(t, numpy.double):
prefix = 'd'
elif issubclass(t, numpy.longdouble):
prefix = 'd'
elif issubclass(t, numpy.csingle):
prefix = 'c'
elif issubclass(t, numpy.cdouble):
prefix = 'z'
elif issubclass(t, numpy.clongdouble):
prefix = 'z'
else:
prefix = 'd'
return prefix
prefix_to_order = dict(s=3, d=2, c=1, z=0)
order_to_prefix = ['s', 'd', 'c', 'z']
prefix_to_dtype = dict(s=numpy.single, d=numpy.double,
c=numpy.csingle, z=numpy.cdouble)
def find_best_lapack_type(arrays):
if not arrays:
return 'd', numpy.double, False
ordering = []
for i in range(len(arrays)):
t = arrays[i].dtype.type
prefix = cast_to_lapack_prefix(t)
order = prefix_to_order[prefix]
ordering.append((order, prefix, i))
ordering.sort()
_, required_prefix, lowest_array_index = ordering[0]
dtype = prefix_to_dtype[required_prefix]
isfortran = numpy.isfortran(arrays[lowest_array_index])
return required_prefix, dtype, isfortran
def get_lapack_funcs(names, arrays=()):
"""Return available LAPACK function objects with names.
arrays are used to determine the optimal prefix of
LAPACK routines.
"""
#If force_clapack is True then available Atlas routine
#is returned for column major storaged arrays with
#rowmajor argument set to False.
force_clapack=False #XXX: Don't set it true! The feature is unreliable
# and may cause incorrect results.
# See test_basic.test_solve.check_20Feb04_bug.
required_prefix, dtype, isfortran = find_best_lapack_type(arrays)
# Default lookup:
if isfortran:
# prefer Fortran code for leading array with column major order
m1, m2 = flapack, clapack
else:
# in all other cases, C code is preferred
m1, m2 = clapack, flapack
if not _use_force_clapack:
force_clapack = False
funcs = []
m1_name = m1.__name__.split('.')[-1]
m2_name = m2.__name__.split('.')[-1]
for name in names:
func_name = required_prefix + name
func = getattr(m1,func_name,None)
if func is None:
func = getattr(m2,func_name)
func.module_name = m2_name
else:
func.module_name = m1_name
if force_clapack and m1 is flapack:
func2 = getattr(m2,func_name,None)
if func2 is not None:
exec _colmajor_func_template % {'func_name':func_name}
func = new.function(func_code,
{'clapack_func':func2},
func_name)
func.module_name = m2_name
func.__doc__ = func2.__doc__
func.prefix = required_prefix
func.dtype = dtype
funcs.append(func)
return tuple(funcs)
_colmajor_func_template = '''\
def %(func_name)s(*args,**kws):
if "rowmajor" not in kws:
kws["rowmajor"] = 0
return clapack_func(*args,**kws)
func_code = %(func_name)s.func_code
'''
|
