|
2 | 2 |
|
3 | 3 | import cython |
4 | 4 | from cpython.mem cimport PyMem_Malloc, PyMem_Free |
5 | | -from libc.string cimport memcpy |
6 | | - |
7 | 5 | from scipy.linalg.cython_lapack cimport sgetrf, dgetrf, cgetrf, zgetrf |
8 | 6 | from scipy.linalg._cythonized_array_utils cimport lapack_t, swap_c_and_f_layout |
9 | 7 |
|
@@ -103,7 +101,9 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a, |
103 | 101 | # as final. Solution without argsort : ipiv[perm] = np.arange(m) |
104 | 102 | for ind1 in range(m): |
105 | 103 | ipiv[perm[ind1]] = ind1 |
106 | | - memcpy(&perm[0], ipiv, m*sizeof(int)) |
| 104 | + for ind1 in range(m): |
| 105 | + perm[ind1] = ipiv[ind1] |
| 106 | + |
107 | 107 | finally: |
108 | 108 | PyMem_Free(ipiv) |
109 | 109 |
|
@@ -135,24 +135,22 @@ cdef void lu_decompose(cnp.ndarray[lapack_t, ndim=2] a, |
135 | 135 | # rows from b as dictated by perm |
136 | 136 |
|
137 | 137 | if m > n: |
138 | | - memcpy(bb, &a[0, 0], m*mn*sizeof(lapack_t)) |
| 138 | + b[:, :] = a[:, :] |
| 139 | + # memcpy(bb, &a[0, 0], m*mn*sizeof(lapack_t)) |
139 | 140 | for ind1 in range(m): |
140 | 141 | if perm[ind1] == ind1: |
141 | 142 | continue |
142 | 143 | else: |
143 | | - memcpy(&a[ind1, 0], |
144 | | - bb + (perm[ind1]*mn), # row stride |
145 | | - mn*sizeof(lapack_t)) # copy one row of memory |
| 144 | + a[ind1, :] = b[perm[ind1], :] |
146 | 145 |
|
147 | 146 | else: # same but for lu array |
148 | | - memcpy(bb, &lu[0, 0], mn*n*sizeof(lapack_t)) |
| 147 | + b[:mn, :mn] = lu[:, :] |
| 148 | + # memcpy(bb, &lu[0, 0], mn*n*sizeof(lapack_t)) |
149 | 149 | for ind1 in range(mn): |
150 | 150 | if perm[ind1] == ind1: |
151 | 151 | continue |
152 | 152 | else: |
153 | | - memcpy(&lu[ind1, 0], |
154 | | - bb + (perm[ind1]*mn), |
155 | | - mn*sizeof(lapack_t)) |
| 153 | + lu[ind1, :] = b[perm[ind1], :mn] |
156 | 154 |
|
157 | 155 |
|
158 | 156 | @cython.nonecheck(False) |
|
0 commit comments