[Cython] Non-deterministic behavoiur?

[Dave Hirschfeld]

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]
#