Mat.getDenseArray unsafe

[Lawrence Mitchell]

getDenseArray is implemented as:

    def getDenseArray(self):
        cdef PetscInt m=0, N=0, lda=0
        cdef PetscScalar *data = NULL
        CHKERR( MatGetLocalSize(self.mat, &m, NULL) )
        CHKERR( MatGetSize(self.mat, NULL, &N) )
        lda = m # CHKERR( MatDenseGetLDA(self.mat, &ld) )
        CHKERR( MatDenseGetArray(self.mat, &data) )
        cdef int typenum = NPY_PETSC_SCALAR
        cdef int itemsize = <int>sizeof(PetscScalar)
        cdef int flags = NPY_ARRAY_FARRAY
        cdef npy_intp dims[2], strides[2]
        dims[0] = <npy_intp>m; strides[0] = <npy_intp>sizeof(PetscScalar);
        dims[1] = <npy_intp>N; strides[1] = <npy_intp>(lda*sizeof(PetscScalar));
        array = <object>PyArray_New(<PyTypeObject*>ndarray, 2, dims, typenum,
                                    strides, data, itemsize, flags, NULL)
        CHKERR( MatDenseRestoreArray(self.mat, &data) )
        return array

That is, we get a reference to the array inside the matrix with MatDenseGetArray, create a numpy array wrapping that data, then call MatDenseRestoreArray. But now, there are two different ideas of the scope of the array. If I do:

mat = createDense(...)

arr = mat.getDenseArray()

del mat

Then arr is now pointing to deallocated memory, because it holds no reference to the matrix it came from. I guess to do this "properly" one would have to mimic the _Vec_buffer stuff for Vecs? Alternately, this could just return a copy of the data (as happens in various other places). WDYT?

[Lisandro Dalcin]

Your diagnosis is correct. The Mat.getDenseArray() is unsafe. As using dense matrices in the context of PETSc is a bit unusual (at least from my POV), I did not care too much about the potential issues.

One obvious fix is to make a copy, but given that dense matrices are quite "fat" things, I think this would be bad. So the reasonable alternative is to implement something similar to _Vec_buffer, let say _MatDense_buffer.