Using the following test code: import numpy as np from lapack import dgelsy from numpy.linalg import lstsq A = np.array([[ 0.12, -8.19, 7.69, -2.26, -4.71], [-6.91, 2.22, -5.12, -9.08, 9.96], [-3.33, -8.94, -6.72, -4.4 , -9.98], [ 3.97, 3.33, -2.74, -7.92, -3.2 ]]) # b = np.array([[ 7.3 , 0.47, -6.28], [ 1.33, 6.58, -3.42], [ 2.68, -1.71, 3.46], [-9.62, -0.79, 0.41]]) # print '#################' print '# ACTUAL RESULT #' print '#################' print lstsq(A, b)[0] print print '#################' print '# DGELSY RESULT #' print '#################' print dgelsy(A, b) I get: ################# # ACTUAL RESULT # ################# [[-0.6858 -0.2431 0.0642] [-0.795 -0.0836 0.2118] [ 0.3767 0.1208 -0.6541] [ 0.2885 -0.2415 0.4176] [ 0.2916 0.3525 -0.3015]] ################# # DGELSY RESULT # ################# [[-0.6858 -0.2431 0.0642] [-0.795 -0.0836 0.2118] [ 0.3767 0.1208 -0.6541] [ 0.2885 -0.2415 0.4176] [ 0.2916 0.3525 -0.3015]] All well and good, however if I type the `tmp` variable as a memview in the following code cdef double[:] res cdef double[:,:] tmp tmp = np.zeros([ldb, nrhs], order='F', dtype=np.float64) tmp[:b.shape[0]] = b res = np.ravel(tmp, order='F') the result changes?!? ################# # DGELSY RESULT # ################# [[-0.7137 -0.2429 0.0876] [-0.2882 -0.0884 -0.2117] [-0.4282 0.1284 0.0185] [ 0.9564 -0.2478 -0.1404] [ 0.3625 0.3519 -0.3607]] Remove the `cdef double[:,:] tmp` and I'm back to the correct result. Does this make any sense? To try and figure out what was going on I put in a couple of debugging print statements: print 'res = ', repr(np.asarray(res)) print 'res.flags = {{{flags}}}'.format(flags=np.asarray(res).flags) Only changing these lines resulted in the same incorrect result ################# # DGELSY RESULT # ################# res = array([ 7.3 , 1.33, 2.68, -9.62, 0., 0.47, 6.58, -1.71, -0.79, 0., -6.28, -3.42, 3.46, 0.41, 0.]) res.flags = { C_CONTIGUOUS : True F_CONTIGUOUS : True OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False} [[-0.7137 -0.2429 0.0876] [-0.2882 -0.0884 -0.2117] [-0.4282 0.1284 0.0185] [ 0.9564 -0.2478 -0.1404] [ 0.3625 0.3519 -0.3607]] Removing (only) the print statements again gave me the correct results. So, it seems either typing the array as a memview or printing res will screw up the calculation. The cython code is given below. Any ideas if this is a cython bug or something I'm doing wrong? Thanks, Dave cdef extern from "mkl_lapack.h" nogil: void DGELSY(const MKL_INT* m, const MKL_INT* n, const MKL_INT* nrhs, double* a, const MKL_INT* lda, double* b, const MKL_INT* ldb, MKL_INT* jpvt, const double* rcond, MKL_INT* rank, double* work, const MKL_INT* lwork, MKL_INT* info) @cython.embedsignature(True) @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) def dgelsy(double[:,:] A, double[:,:] b, double rcond=1e-15, overwrite_A=False, overwrite_b=False): cdef MKL_INT rank, info cdef MKL_INT m = A.shape[0] cdef MKL_INT n = A.shape[1] cdef MKL_INT nrhs = b.shape[1] cdef MKL_INT lda = m cdef MKL_INT ldb = max(m, n) cdef MKL_INT lwork = -1 cdef double worksize = 0 #cdef double[:,:] tmp cdef double[:] res, work cdef MKL_INT[:] jpvt if b.shape[0] != m: message = "b.shape[0] must be equal to A.shape[0].\n" message += "b.shape[0] = {b.shape[0]}\n" message += "A.shape[0] = {A.shape[0]}\n" raise MKL_LAPACK_ERROR(message.format(A=A, b=b)) flags = np.asarray(A).flags if not flags['F_CONTIGUOUS'] or not overwrite_A: A = A.copy_fortran() flags = np.asarray(b).flags if not flags['F_CONTIGUOUS'] or not overwrite_b or b.shape[0] < n: tmp = np.zeros([ldb, nrhs], order='F', dtype=np.float64) tmp[:b.shape[0]] = b res = np.ravel(tmp, order='F') else: res = np.ravel(b, order='F') #print 'res = ', repr(np.asarray(res)) #print 'res.flags = {{{flags}}}'.format(flags=np.asarray(res).flags) jpvt = np.empty(n, dtype=np.int32) DGELSY(&m, &n, &nrhs, &A[0,0], &lda, &res[0], &ldb, &jpvt[0], &rcond, &rank, &worksize, &lwork, &info) if info != 0: message = "Parameter {i} had an illegal value when calling gelsy." raise MKL_LAPACK_ERROR(message.format(i=info)) lwork = int(worksize) work = np.empty(lwork, dtype=np.float64) DGELSY(&m, &n, &nrhs, &A[0,0], &lda, &res[0], &ldb, &jpvt[0], &rcond, &rank, &worksize, &lwork, &info) if info != 0: message = "Parameter {i} had an illegal value when calling gelsy." raise MKL_LAPACK_ERROR(message.format(i=info)) return np.asarray(res).reshape(nrhs, -1).T[:n] #