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]
#
_______________________________________________
cython-devel mailing list
cython-devel@python.org
http://mail.python.org/mailman/listinfo/cython-devel
