How to identify memory leaks in extensions?

[Tom Loredo]

Hi folks-

One of my extensions using NumPy arrays appears to have a memory leak.
Do any of you more experience extension programs have any tips and
techniques you'd care to share on identifying the source of memory
leaks in a dynamically loaded module?  I'm on a Mac and Solaris, by
the way.  I suspect the problem is that I haven't DECREF'd an 
ArrayObject correctly, but I don't find the NumPy documentation
very helpful with advice on when to DECREF and INCREF (e.g., the
"trace" and "matrix_vector" examples in Ch. 12 appear to my
nonexpert eye to offer inconsistent advice on whether to DECREF
an array created by PyArray_PyContiguousFromObject and only
used within the extension).

For the curious, the code is posted below if you'd like to offer
specific rather than general advice.  Either is appreciated!

Thanks in advance for any help you can offer,
Tom Loredo

/*------------------------------------------------------------------------------
* Calculate the model metric & related quantities needed for the marginal
* distribution for the nonlinear parameters.  Return them in a tuple:
*
* (metric, L, rdet, P, A, S)
*
*    metric = model metric (MxM matrix)
*    L      = Cholesky factorization of metric (lower triangular MxM)
*    rdet   = Sqrt of determinant of the metric
*    P      = Projections of data on models (M vector)
*    A      = Amplitude estimates (M vector)
*    S      = Sufficient statistic
*
* The Cholesky algorithms are adapted from *Numerical Recipes*.
------------------------------------------------------------------------------*/
static char vba_metricAnalysis_doc[] =
"metricAnalysis(dmat,data):  Calculate the metric from the design matrix; get
various derived quantities from it & the data via Cholesky.";

static PyObject *
vba_metricAnalysis (PyObject *self, PyObject *args) {
    int mdims[2], pdims[1];
    int nmod, nsamp, a, b, i, m;
    PyObject *input1, *input2;
    PyArrayObject *dmat, *smpls;
    PyArrayObject *metric, *cfmetric, *proj, *amps;
    double *dmdata, *mdata, *cfdata, *diag, *pdata, *sdata, *adata;
    double sum, rdet;

/** Parse the design matrix & data; make them contiguous. */
    Py_TRY(PyArg_ParseTuple(args, "OO", &input1, &input2));
    dmat = (PyArrayObject *)
            PyArray_ContiguousFromObject(input1, PyArray_DOUBLE, 2, 2);
    if (dmat == NULL) return NULL;
    smpls = (PyArrayObject *)
            PyArray_ContiguousFromObject(input2, PyArray_DOUBLE, 1, 1);
    if (smpls == NULL) return NULL;

/** Check that their dimensions are consistent. */
    if (dmat->dimensions[1] != smpls->dimensions[0]) {
        PyErr_SetString(PyExc_ValueError, "inconsistent data dimensions!");
            return NULL;
    }

/** Store the # of models & samples; make vector references to dmat & data. */
    nmod = dmat->dimensions[0];
    nsamp = dmat->dimensions[1];
    dmdata = DDATA(dmat);
    sdata = DDATA(smpls);

/** Calculate the metric. */
    mdims[0] = nmod;
    mdims[1] = nmod;
    metric = (PyArrayObject *)PyArray_FromDims(2, mdims, PyArray_DOUBLE);
    if (metric == NULL) return NULL;
    mdata = DDATA(metric);
    for (a=0; a<nmod; a++) {
        for (b=0; b<nmod; b++) {
            sum = 0.;
            for (i=0; i<nsamp; i++) {
                sum += *(dmdata+a*nsamp+i) * *(dmdata+b*nsamp+i);
            }
            *(mdata+a*nmod+b) = sum;
        }
    }

/** Calculate the projections of the data on the models. */
    pdims[0] = nmod;
    proj = (PyArrayObject *)PyArray_FromDims(1, pdims, PyArray_DOUBLE);
    pdata = DDATA(proj);
    for (a=0; a<nmod; a++) {
        sum = 0.;
        for (i=0; i<nsamp; i++) {
            sum += *(dmdata+a*nsamp+i) * *(sdata+i);
        }
        *(pdata+a) = sum;
    }

/** Perform a Cholesky factorization of the metric. Calculate
*** the sqrt(determinant) along the way.  */
    cfmetric = (PyArrayObject *)PyArray_Copy(metric);
    if (cfmetric == NULL) return NULL;
    cfdata = DDATA(cfmetric);
    diag = (double *) malloc(nmod*sizeof(double));
    rdet = 1.;
    for (a=0; a<nmod; a++) {
        for (b=a; b<nmod; b++) {
            sum = *(cfdata+a*nmod+b);
            for (m=a-1; m>=0; m--) {
                sum -= *(cfdata+a*nmod+m) * *(cfdata+b*nmod+m);
            }
            if (a == b) {
                if (sum <= 0.) {
                    PyErr_SetString(PyExc_ValueError,
                        "metric is not positive definite!");
                    return NULL;
                }
                *(diag+a) = sqrt(sum);
                rdet = rdet * *(diag+a);
            } else *(cfdata+b*nmod+a) = sum/ *(diag+a);
        }
    }

/** We have the factorization in the lower triangle and diag;
*** replace the diagonal and zero the upper triangle. */
    for (a=0; a<nmod; a++) {
        for (b=a; b<nmod; b++) {
            if (a == b) {
                *(cfdata+a*nmod+a) = *(diag+a);
            } else *(cfdata+a*nmod+b) = 0.;
        }
    }

/** Calculate the amplitude estimates by backsubstitution. */
    pdims[0] = nmod;
    amps = (PyArrayObject *)PyArray_FromDims(1, pdims, PyArray_DOUBLE);
    adata = DDATA(amps);
    for (a=0; a<nmod; a++) {
        sum = *(pdata+a);
        for (b=a-1; b>=0; b--) {
            sum -= *(cfdata+a*nmod+b) * *(adata+b);
        }
        *(adata+a) = sum / *(cfdata+a*nmod+a);
    }
    for (a=nmod-1; a>=0; a--) {
        sum = *(adata+a);
        for (b=a+1; b<nmod; b++) {
            sum -= *(cfdata+b*nmod+a) * *(adata+b);
        }
        *(adata+a) = sum / *(cfdata+a*nmod+a);
    }

/** Calculate the sufficient statistic. */
    sum = 0.;
    for (a=0; a<nmod; a++) {
        sum += *(adata+a) * *(pdata+a);
    }

/** Clean up and return everything in a tuple. */
    free(diag);
    Py_DECREF(dmat);
    Py_DECREF(smpls);
    return Py_BuildValue("OOfOOf",metric,cfmetric,rdet,proj,amps,sum);
}

/*------------------------------------------------------------------------------
* Calculate covariance matrix for the amplitudes from the Cholesky
* factorization of the metric.
------------------------------------------------------------------------------*/
static char vba_covar_doc[] =
"covar(L):  Calculate the covariance matrix for the amplitudes from the Cholesky
factorization of the metric.";

static PyObject *
vba_covar (PyObject *self, PyObject *args) {
    int nmod, a, b, m;
    PyObject *input;
    PyArrayObject *cfmetric, *covar;
    double *cfdata, *cvdata;
    double sum;

/** Parse the Cholesky factor & make it contiguous. */
    Py_TRY(PyArg_ParseTuple(args, "O", &input));
    cfmetric = (PyArrayObject *)
            PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2);
    if (cfmetric == NULL) return NULL;

/** Copy the matrix and gather needed info. */
    covar = (PyArrayObject *)PyArray_Copy(cfmetric);
    if (covar == NULL) return NULL;
    nmod = cfmetric->dimensions[0];
    cvdata = DDATA(covar);

/** Find the inverse of the metric to use as covariances.  First find
*** the inverse of L; then multiply in place for the full upper inverse;
*** then fill in the lower triangle. */
    for (a=0; a<nmod; a++) {
        *(cvdata+a*nmod+a) = 1. / *(cvdata+a*nmod+a);
        for (b=a+1; b<nmod; b++) {
            sum = 0.;
            for (m=a; m<b; m++) sum -= *(cvdata+b*nmod+m) * *(cvdata+m*nmod+a);
            *(cvdata+b*nmod+a) = sum / *(cvdata+b*nmod+b);
        }
    }
    for (a=0; a<nmod; a++) {
        for (b=a; b<nmod; b++) {
            sum = 0.;
            for (m=b; m<nmod; m++) sum += *(cvdata+m*nmod+a) *
*(cvdata+m*nmod+b);
            *(cvdata+a*nmod+b) = sum;
        }
    }
    for (a=0; a<nmod; a++) {
        for (b=a+1; b<nmod; b++) {
            *(cvdata+b*nmod+a) = *(cvdata+a*nmod+b);
        }
    }

/** Clean up and return the covariance matrix. */
    Py_DECREF(cfmetric);
    return PyArray_Return(covar);
}