Hi list,

I'm suffering from a strange segfault and would appreciate your help.

I'm calling a small C function using ctypes / numpy.ctypeslib. The function 
works in the sense that it returns correct results. After calling the 
function however I can reliably evoke a segfault by using readline tab 
completion.

I'm not very experienced, but this smells like a memory management bug to me, 
which is strange, because I'm not doing any mallocing/freeing at all in the C 
code.

I could not reproduce the bug in a debug build of python (--without-pymalloc) 
or on another machine. The crashing machine is an eight-way opteron.

Maybe I should mention that the C function calls two lapack fortran functions. 
Can this cause problems?

Anyway, I'm at a loss. Please help!

I've attached the files for reference.

Thanks Martin

P.S.: Here is what valgrind finds in the debug build:

-bash-3.1$ valgrind ~/local/debug/bin/python
==16266== Memcheck, a memory error detector.
==16266== Copyright (C) 2002-2006, and GNU GPL'd, by Julian Seward et al.
==16266== Using LibVEX rev 1658, a library for dynamic binary translation.
==16266== Copyright (C) 2004-2006, and GNU GPL'd, by OpenWorks LLP.
==16266== Using valgrind-3.2.1, a dynamic binary instrumentation framework.
==16266== Copyright (C) 2000-2006, and GNU GPL'd, by Julian Seward et al.
==16266== For more details, rerun with: -v
==16266==
Python 2.5.1 (r251:54863, Aug 24 2007, 16:13:26)
[GCC 4.1.1 20070105 (Red Hat 4.1.1-51)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> execfile ('recttest.py')
--16266-- DWARF2 CFI reader: unhandled CFI instruction 0:10
--16266-- DWARF2 CFI reader: unhandled CFI instruction 0:10
78.6006829739 4
78.6006829739
[92353 refs]
>>>
==16266== Conditional jump or move depends on uninitialised value(s)
==16266==    at 0x41361F: parsetok (parsetok.c:189)
==16266==    by 0x4131B6: PyParser_ParseFileFlags (parsetok.c:89)
==16266==    by 0x4E01D2: PyParser_ASTFromFile (pythonrun.c:1381)
==16266==    by 0x4DE15C: PyRun_InteractiveOneFlags (pythonrun.c:770)
==16266==    by 0x4DDF15: PyRun_InteractiveLoopFlags (pythonrun.c:721)
==16266==    by 0x4DDD70: PyRun_AnyFileExFlags (pythonrun.c:690)
==16266==    by 0x412E05: Py_Main (main.c:523)
==16266==    by 0x411D62: main (python.c:23)
[92353 refs]
[37870 refs]
==16266==
==16266== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 98 from 1)
==16266== malloc/free: in use at exit: 2,791,249 bytes in 17,090 blocks.
==16266== malloc/free: 174,713 allocs, 157,623 frees, 376,495,377 bytes 
allocated.
==16266== For counts of detected errors, rerun with: -v
==16266== searching for pointers to 17,090 not-freed blocks.
==16266== checked 5,156,624 bytes.
==16266==
==16266== LEAK SUMMARY:
==16266==    definitely lost: 0 bytes in 0 blocks.
==16266==      possibly lost: 35,400 bytes in 103 blocks.
==16266==    still reachable: 2,755,849 bytes in 16,987 blocks.
==16266==         suppressed: 0 bytes in 0 blocks.
==16266== Reachable blocks (those to which a pointer was found) are not shown.
==16266== To see them, rerun with: --show-reachable=yes


#include <stdlib.h>
int dgetrf_ (int *, int *, double *, int *, int *, int *);
void dgetrs_ (char *, int *, int *, double *, int *, int *, double *, int *,
	      int *);

#define RECTIFIED_SOLVER_MAX_ITER -1
#define RECTIFIED_SOLVER_FACTORISE -2
#define RECTIFIED_SOLVER_SOLVE -3

/* Copy full data to W, LU, s and iWs. Initialise inds to 0, 1, 2, 3, ...
   Leaks inds_rp for potential use with MIXED_QUADRANT. */
#define INIT								\
{									\
    for (i = 0, inds_rp = inds, LU_p = LU, W_p = W, W_rp = W_orig,	\
	     iWs_p = iWs, s_p = s, s_rp = s_orig;			\
	 i < n; i++, inds_rp++, iWs_p++, s_p++, s_rp++)			\
    {									\
	*inds_rp = i;							\
	*iWs_p = *s_p = *s_rp;						\
	for (next = W_rp + n; W_rp < next; LU_p++, W_p++, W_rp++)	\
	    *LU_p = *W_p = *W_rp;					\
    }									\
}

/* Apply mask *inplace* to weight matrix W, stimulus s and index table inds.
   Overwrite n with new size.
   Simultaneously copy W to LU and s to iWs.
   LAPACK will overwrite these copies.
   Leaks inds_rp for potential use with MIXED_QUADRANT. */
#define APPLY_REL_MASK							  \
{									  \
    for (y_mask_p = mask, s_p = s_rp = s, iWs_p = iWs,			  \
	     inds_p = inds_rp = inds, W_p = W_rp = W, LU_p = LU;	  \
	 y_mask_p < mask_top;						  \
	 y_mask_p++, s_p++, inds_p++)					  \
	if (*y_mask_p)							  \
	{								  \
	    *iWs_p = *s_rp = *s_p;					  \
	    iWs_p++; s_rp++;						  \
	    *inds_rp = *inds_p;						  \
	    inds_rp++;							  \
	    for (x_mask_p = mask; x_mask_p < mask_top; x_mask_p++, W_p++) \
		if (*x_mask_p)						  \
		{							  \
		    *LU_p = *W_rp = *W_p;				  \
		    LU_p++; W_rp++;					  \
		}							  \
	}								  \
	else								  \
	    W_p += n;							  \
    n = s_rp - s;							  \
}

/* Copy temp_inds to inds and apply them
   to full data to get weight matrix W, stimulus s.
   Simultaneously copy W to LU and s to iWs.
   LAPACK will overwrite these copies.
   Uses inds_p leaked from MIXED_QUADRANT. */
#define APPLY_ABS_MASK					\
{							\
    n = inds_p - temp_inds;				\
    for (inds_rp = inds, y_inds_p = temp_inds,		\
	     LU_p = LU, W_p = W, iWs_p = iWs, s_p = s;	\
	 y_inds_p < inds_p;				\
	 inds_rp++, y_inds_p++, iWs_p++, s_p++)		\
    {							\
	*inds_rp = *y_inds_p;				\
	*iWs_p = *s_p = s_orig [*y_inds_p];		\
	W_rp = W_orig + *y_inds_p * n_orig;		\
	for (x_inds_p = temp_inds; x_inds_p < inds_p;	\
	     x_inds_p++, LU_p++, W_p++)			\
	    *LU_p = *W_p = W_rp [*x_inds_p];		\
    }							\
}

/* Uses inds_rp leaked from APPLY_REL_MASK
   and leaks inds_p for APPLY_ABS_MASK. */
#define MIXED_QUADRANT							\
{									\
    done = 1;								\
    for (y_inds_p = inds, inds_p = temp_inds, j = 0; j < n_orig; j++)	\
	if (j == *y_inds_p)						\
	{								\
	    *inds_p = j;						\
	    inds_p++; y_inds_p++;					\
	}								\
	else								\
	{								\
	    W_p = W_orig + j * n_orig;					\
	    for (sum = 0.0, x_inds_p = inds, iWs_p = iWs;		\
		 x_inds_p < inds_rp; x_inds_p++, iWs_p++)		\
		sum += W_p [*x_inds_p] * *iWs_p;			\
	    if (sum < s_orig [j])					\
	    {								\
		done = 0;						\
		*inds_p = j;						\
		inds_p++;						\
	    }								\
	}								\
}

#define SOLVE							\
{								\
    dgetrf_ (&n, &n, LU, &n, perm, info);			\
    if (*info != 0)						\
	return RECTIFIED_SOLVER_FACTORISE;			\
    dgetrs_ ("T", &n, &p1, LU, &n, perm, iWs, &n, info);	\
    if (*info != 0)						\
	return RECTIFIED_SOLVER_SOLVE;				\
}

#define CHECK_SOL				\
{						\
    mask_top = mask + n;			\
    done = 1;					\
    for (x_mask_p = mask, iWs_p = iWs;		\
	 x_mask_p < mask_top;			\
         x_mask_p++, iWs_p++)			\
	if (*iWs_p > 0.0)			\
	    *x_mask_p = 1;			\
	else					\
	{					\
	    done = 0;				\
	    *x_mask_p = 0;			\
	}					\
}

#define SUCCESS					\
{						\
    *info = i;					\
    return n;					\
}

/* W, LU, s, iWs, inds, temp_inds, mask, perm must be allocated by user.
   On success will return the number of non silent cells.
   W will contain the reduced matrix, LU and perm the pivoted LU factorisation,
   s the reduced stimulus, iWs the reduced solution vector
   and inds the indices of non silent cells.
   Upon failure a negative error code is returned,
   if the error comes from LAPACK the LAPACK error code is in perm [n_orig] */
int solve (int n_orig, double *W_orig, double *s_orig, double *W, double *LU,
	   double *s, double *iWs, ssize_t *inds, ssize_t *temp_inds,
	   char *mask, int *perm, int max_iter)
{
    ssize_t *inds_p, *inds_rp, *x_inds_p, *y_inds_p;
    int i, j, n = n_orig, *info = perm + n, p1 = 1;
    char done, *x_mask_p, *y_mask_p, *mask_top;
    double sum, *s_p, *W_p, *s_rp, *W_rp, *iWs_p, *LU_p, *next;

    INIT;
    for (i = max_iter; i > 0; i--)
    {
	while (1)
	{
	    SOLVE;
	    CHECK_SOL;
	    if (done)
		break;
	    APPLY_REL_MASK;
	}
	MIXED_QUADRANT;
	if (done)
	    SUCCESS;
	APPLY_ABS_MASK;
    }
    return RECTIFIED_SOLVER_MAX_ITER;
}

Attachment: cMontecarlo.py
Description: application/python

_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Reply via email to