diff options
Diffstat (limited to 'system/type/pdp_matrix.c')
| -rw-r--r-- | system/type/pdp_matrix.c | 649 |
1 files changed, 0 insertions, 649 deletions
diff --git a/system/type/pdp_matrix.c b/system/type/pdp_matrix.c deleted file mode 100644 index 1334864..0000000 --- a/system/type/pdp_matrix.c +++ /dev/null @@ -1,649 +0,0 @@ - -#include "pdp.h" -#include "pdp_internals.h" - -/* the class object */ -static t_pdp_class *matrix_class; - - - -/* declare header and subheader variables and exit with errval if invalid */ -#define VALID_MATRIX_HEADER(packet, header, subheader, mtype, errval) \ -t_pdp * header = pdp_packet_header( packet ); \ -t_matrix * subheader = (t_matrix *)pdp_packet_subheader( packet ); \ -if (! header ) return errval; \ -if (PDP_MATRIX != header->type) return errval; \ -if (mtype) {if (subheader->type != mtype) return errval;} - - -int pdp_packet_matrix_isvalid(int p) -{ - VALID_MATRIX_HEADER(p, h, m, 0, 0); - return 1; -} - - -void *pdp_packet_matrix_get_gsl_matrix(int p, u32 type) -{ - VALID_MATRIX_HEADER(p, h, m, type, 0); - return &m->matrix; -} - -void *pdp_packet_matrix_get_gsl_vector(int p, u32 type) -{ - VALID_MATRIX_HEADER(p, h, m, type, 0); - return &m->vector; -} - -int pdp_packet_matrix_get_type(int p) -{ - VALID_MATRIX_HEADER(p, h, m, 0, 0); - return m->type; -} - -int pdp_packet_matrix_isvector(int p) -{ - VALID_MATRIX_HEADER(p, h, m, 0, 0); - return ((m->rows == 1) || (m->columns == 1)); -} - -int pdp_packet_matrix_ismatrix(int p) -{ - VALID_MATRIX_HEADER(p, h, m, 0, 0); - return ((m->rows != 1) && (m->columns != 1)); -} - - - - -/* gsl blas matrix/vector multiplication: - -vector.vector - - (const gsl_vector * x, const gsl_vector * y, double * result) - -gsl_blas_sdot -gsl_blas_ddot -gsl_blas_cdot -gsl_blas_zdot -gsl_blas_cdotu -gsl_blas_zdotu - -matrix.vector - ( - CBLAS_TRANSPOSE_t - TransA, - double alpha, - const gsl_matrix * A, - const gsl_vector * x, - double beta, - gsl_vector * y - ) - -gsl_blas_sgemv -gsl_blas_dgemv -gsl_blas_cgemv -gsl_blas_zgemv - -matrix.matrix - ( - CBLAS_TRANSPOSE_t TransA, - CBLAS_TRANSPOSE_t TransB, - double alpha, - const gsl_matrix * A, - const gsl_matrix * B, - double beta, - gsl_matrix * C - ) - -gsl_blas_sgemm -gsl_blas_dgemm -gsl_blas_cgemm -gsl_blas_zgemm - -*/ - -/* compute the matrix inverse using the LU decomposition */ -/* it only works for double real/complex */ -int pdp_packet_matrix_LU_to_inverse(int p) -{ - int new_p; - u32 n; - int type = pdp_packet_matrix_get_type(p); - t_matrix *sheader = (t_matrix *)pdp_packet_subheader(p); - gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type); - gsl_matrix *im; - if (!m) return -1; - n = m->gsl_rows; - if (n != m->gsl_columns) return -1; - new_p = pdp_packet_new_matrix(n, n, type); - if (-1 == new_p) return -1; - im = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(new_p, type); - - switch(type){ - case PDP_MATRIX_TYPE_RDOUBLE: - gsl_linalg_LU_invert (m, &sheader->perm, im); break; - case PDP_MATRIX_TYPE_CDOUBLE: - gsl_linalg_complex_LU_invert ((gsl_matrix_complex *)m, &sheader->perm, (gsl_matrix_complex *)im); break; - default: - pdp_packet_mark_unused(new_p); - new_p = -1; - } - return new_p; -} -/* compute the LU decomposition of a square matrix */ -/* it only works for double real/complex */ -int pdp_packet_matrix_LU(int p) -{ - int p_LU, bytes; - u32 n; - int type = pdp_packet_matrix_get_type(p); - t_matrix *sh_m_LU; - t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p); - gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type); - gsl_matrix *m_LU; - if (!m) return -1; - n = m->gsl_rows; - if (n != m->gsl_columns) return -1; - p_LU = pdp_packet_new_matrix(n, n, type); - if (-1 == p_LU) return -1; - sh_m_LU = (t_matrix *)pdp_packet_subheader(p_LU); - m_LU = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_LU, type); - /* copy matrix data: move this to copy method */ - memcpy(pdp_packet_data(p_LU), pdp_packet_data(p), sh_m->block.size); - - switch(type){ - case PDP_MATRIX_TYPE_RDOUBLE: - gsl_linalg_LU_decomp (m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break; - case PDP_MATRIX_TYPE_CDOUBLE: - gsl_linalg_complex_LU_decomp ((gsl_matrix_complex *)m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break; - default: - pdp_packet_mark_unused(p_LU); - p_LU = -1; - } - return p_LU; -} - -int pdp_packet_matrix_LU_solve(int p_matrix, int p_vector) -{ - int type = pdp_packet_matrix_get_type(p_matrix); - int p_result_vector = pdp_packet_clone_rw(p_vector); - gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_matrix, type); - gsl_vector *v = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_vector, type); - gsl_vector *rv = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_result_vector, type); - t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p_matrix); - - if (!(m && v && rv)) goto error; - - switch(type){ - case PDP_MATRIX_TYPE_RDOUBLE: - if (gsl_linalg_LU_solve (m, &sh_m->perm, v, rv)) goto error; break; - case PDP_MATRIX_TYPE_CDOUBLE: - if(gsl_linalg_complex_LU_solve ((gsl_matrix_complex*)m, &sh_m->perm, - (gsl_vector_complex *)v, (gsl_vector_complex *)rv)) goto error; break; - default: - goto error; - } - return p_result_vector; - - error: - pdp_packet_mark_unused(p_result_vector); - //post("error"); - return -1; -} - -/* matrix matrix mul: C is defining type - returns 0 on success */ -int pdp_packet_matrix_blas_mm -( - CBLAS_TRANSPOSE_t TransA, - CBLAS_TRANSPOSE_t TransB, - int pA, - int pB, - int pC, - float scale_r, - float scale_i -) -{ - gsl_complex_float cf_scale = {{scale_r, scale_i}}; - gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}}; - gsl_complex_float cf_one = {{1.0f, 0.0f}}; - gsl_complex cd_one = {{1.0, 0.0}}; - gsl_matrix *mA, *mB, *mC; - int type; - type = pdp_packet_matrix_get_type(pC); - mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type); - mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB,type); - mC = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pC,type); - - if (!(mA && mB)) return 1; - - - switch(type){ - case PDP_MATRIX_TYPE_RFLOAT: - return gsl_blas_sgemm(TransA, TransB, scale_r, (gsl_matrix_float *)mA, - (gsl_matrix_float *)mB, 1.0f, (gsl_matrix_float *)mC); - case PDP_MATRIX_TYPE_RDOUBLE: - return gsl_blas_dgemm(TransA, TransB, (double)scale_r, (gsl_matrix *)mA, - (gsl_matrix *)mB, 1.0, (gsl_matrix *)mC); - case PDP_MATRIX_TYPE_CFLOAT: - return gsl_blas_cgemm(TransA, TransB, cf_scale, (gsl_matrix_complex_float *)mA, - (gsl_matrix_complex_float *)mB, cf_one, (gsl_matrix_complex_float *)mC); - case PDP_MATRIX_TYPE_CDOUBLE: - return gsl_blas_zgemm(TransA, TransB, cd_scale, (gsl_matrix_complex *)mA, - (gsl_matrix_complex *)mB, cd_one, (gsl_matrix_complex *)mC); - default: - return 0; - } -} - -/* matrix vector mul: C is defining type - returns 0 on success */ -int pdp_packet_matrix_blas_mv -( - CBLAS_TRANSPOSE_t TransA, - int pA, - int pb, - int pc, - float scale_r, - float scale_i -) -{ - gsl_complex_float cf_scale = {{scale_r, scale_i}}; - gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}}; - gsl_complex_float cf_one = {{1.0f, 0.0f}}; - gsl_complex cd_one = {{1.0, 0.0}}; - gsl_matrix *mA; - gsl_vector *vb, *vc; - int type; - type = pdp_packet_matrix_get_type(pA); - mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type); - vb = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pb,type); - vc = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pc,type); - - if (!(vb && vc)) return 1; - - - switch(type){ - case PDP_MATRIX_TYPE_RFLOAT: - return gsl_blas_sgemv(TransA, scale_r, (gsl_matrix_float *)mA, - (gsl_vector_float *)vb, 1.0f, (gsl_vector_float *)vc); - case PDP_MATRIX_TYPE_RDOUBLE: - return gsl_blas_dgemv(TransA, (double)scale_r, (gsl_matrix *)mA, - (gsl_vector *)vb, 1.0, (gsl_vector *)vc); - case PDP_MATRIX_TYPE_CFLOAT: - return gsl_blas_cgemv(TransA, cf_scale, (gsl_matrix_complex_float *)mA, - (gsl_vector_complex_float *)vb, cf_one, (gsl_vector_complex_float *)vc); - case PDP_MATRIX_TYPE_CDOUBLE: - return gsl_blas_zgemv(TransA, cd_scale, (gsl_matrix_complex *)mA, - (gsl_vector_complex *)vb, cd_one, (gsl_vector_complex *)vc); - default: - return 0; - } -} - - - -t_pdp_symbol *_pdp_matrix_get_description(t_pdp *header) -{ - t_matrix *m = (t_matrix *)(&header->info.raw); - char description[100]; - char *c = description; - int encoding; - - if (!header) return pdp_gensym("invalid"); - else if (!header->desc){ - /* if description is not defined, try to reconstruct it (for backwards compat) */ - if (header->type == PDP_MATRIX){ - - c += sprintf(c, "matrix"); - - switch(m->type){ - case PDP_MATRIX_TYPE_RFLOAT: c += sprintf(c, "/float/real"); break; - case PDP_MATRIX_TYPE_CFLOAT: c += sprintf(c, "/float/complex"); break; - case PDP_MATRIX_TYPE_RDOUBLE: c += sprintf(c, "/double/real"); break; - case PDP_MATRIX_TYPE_CDOUBLE: c += sprintf(c, "/double/complex"); break; - default: - c += sprintf(c, "/unknown"); goto exit; - } - - c += sprintf(c, "/%dx%d", (int)m->rows, (int)m->columns); - - exit: - return pdp_gensym(description); - } - else return pdp_gensym("unknown"); - } - else return header->desc; -} - - - -static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src); -static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src); -static void _pdp_matrix_reinit(t_pdp *dst); -static void _pdp_matrix_cleanup(t_pdp *dst); - -static size_t _pdp_matrix_mdata_byte_size(u32 rows, u32 columns, u32 type) -{ - size_t dsize = rows * columns; - switch (type){ - case PDP_MATRIX_TYPE_RFLOAT: dsize *= sizeof(float); break; - case PDP_MATRIX_TYPE_CFLOAT: dsize *= 2*sizeof(float); break; - case PDP_MATRIX_TYPE_RDOUBLE: dsize *= sizeof(double); break; - case PDP_MATRIX_TYPE_CDOUBLE: dsize *= 2*sizeof(double); break; - default: dsize = 0; - } - return dsize; -} - - -static size_t _pdp_matrix_pdata_vector_size(u32 rows, u32 columns) -{ - return (rows>columns ? rows : columns); -} - - -static void _pdp_matrix_init(t_pdp *dst, u32 rows, u32 columns, u32 type) -{ - int i; - t_matrix *m = (t_matrix *)(&dst->info.raw); - void *d = ((void *)dst) + PDP_HEADER_SIZE; - int matrix_bytesize = _pdp_matrix_mdata_byte_size(rows, columns, type); - int permsize = _pdp_matrix_pdata_vector_size(rows, columns); - - /* set meta data */ - m->type = type; - m->rows = rows; - m->columns = columns; - - /* set the block data */ - m->block.size = matrix_bytesize; - m->block.data = (double *)d; - - /* set the vector view */ - m->vector.size = (1==columns) ? rows : columns; - m->vector.stride = 1; - m->vector.data = (double *)d; - m->vector.block = &m->block; - m->vector.owner = 0; - - /* set the matrix view */ - m->matrix.gsl_rows = rows; - m->matrix.gsl_columns = columns; - m->matrix.tda = columns; - m->matrix.data = (double *)d; - m->matrix.block = &m->block; - m->matrix.owner = 0; - - /* set the permutation object & init */ - m->perm.size = permsize; - m->perm.data = (size_t *)(d + matrix_bytesize); - for(i=0; i<permsize; i++) m->perm.data[i] = i; - m->signum = -1; - - /* init packet header */ - dst->theclass = matrix_class; - dst->desc = _pdp_matrix_get_description(dst); - -} - -static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src) -{ - t_matrix *m = (t_matrix *)(&src->info.raw); - _pdp_matrix_init(dst, m->rows, m->columns, m->type); - -} -static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src) -{ - _pdp_matrix_clone(dst, src); - memcpy(dst + PDP_HEADER_SIZE, src + PDP_HEADER_SIZE, src->size - PDP_HEADER_SIZE); - //post("matrix copy successful"); -} -static void _pdp_matrix_reinit(t_pdp *dst) -{ - /* nothing to do, assuming data is correct */ -} -static void _pdp_matrix_cleanup(t_pdp *dst) -{ - /* no extra memory to free */ -} - -int pdp_packet_new_matrix(u32 rows, u32 columns, u32 type) -{ - t_pdp *header; - int dsize, p; - - /* compute the blocksize for the matrix data */ - /* if 0, something went wrong -> return invalid packet */ - if (!(dsize = _pdp_matrix_mdata_byte_size(rows, columns, type))) return -1; - dsize += sizeof(size_t) * _pdp_matrix_pdata_vector_size(rows, columns); - - p = pdp_packet_new(PDP_MATRIX, dsize); - if (-1 == p) return -1; - header = pdp_packet_header(p); - - _pdp_matrix_init(header, rows, columns, type); - - return p; -} - -int pdp_packet_new_matrix_product_result(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, int pA, int pB) -{ - gsl_matrix *mA, *mB; - u32 type = pdp_packet_matrix_get_type(pA); - - u32 colA, colB, rowA, rowB; - - /* check if A is a matrix */ - if (!type) return -1; - - mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA, type); - mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB, type); - - /* check if A and B are same type */ - if (!mB) return -1; - - /* get dims A */ - if (TransA == CblasNoTrans){ - rowA = mA->gsl_rows; - colA = mA->gsl_columns; - } - else { - rowA = mA->gsl_columns; - colA = mA->gsl_rows; - } - - /* get dims B */ - if (TransB == CblasNoTrans){ - rowB = mB->gsl_rows; - colB = mB->gsl_columns; - } - else { - rowB = mB->gsl_columns; - colB = mB->gsl_rows; - } - - /* check if sizes are compatible */ - if (colA != rowB) return -1; - - /* create new packet */ - return pdp_packet_new_matrix(rowA, colB, type); -} - -void pdp_packet_matrix_setzero(int p) -{ - t_matrix *m; - if (!pdp_packet_matrix_isvalid(p)) return; - m = (t_matrix *) pdp_packet_subheader(p); - memset(m->block.data, 0, m->block.size); -} - - - -/* type conversion programs */ -static int _pdp_packet_convert_matrix_to_greyimage(int packet, t_pdp_symbol *dest_template) -{ - t_pdp *header = pdp_packet_header(packet); - t_matrix *matrix = (t_matrix *)pdp_packet_subheader(packet); - void *data = pdp_packet_data(packet); - s16 *new_data; - u32 c,r, nbelements; - int new_p; - u32 i; - - c = matrix->columns; - r = matrix->rows; - nbelements = c*r; - - new_p = pdp_packet_new_image_grey(c,r); - if (-1 == new_p) return -1; - new_data = (s16 *)pdp_packet_data(new_p); - - /* convert first channel */ - switch(matrix->type){ - case PDP_MATRIX_TYPE_RFLOAT: - for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i] * (float)0x8000); break; - case PDP_MATRIX_TYPE_CFLOAT: //only copy real channel - for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i<<1] * (float)0x8000); break; - case PDP_MATRIX_TYPE_RDOUBLE: - for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i] * (float)0x8000); break; - case PDP_MATRIX_TYPE_CDOUBLE: //only copy real channel - for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i<<1] * (float)0x8000); break; - default: - pdp_packet_mark_unused(new_p); - new_p = -1; - } - return new_p; -} - -static int _pdp_packet_convert_image_to_matrix(int packet, t_pdp_symbol *dest_template, int type) -{ - t_pdp *header = pdp_packet_header(packet); - t_image *image = pdp_packet_image_info(packet); - s16 *data = (s16 *)pdp_packet_data(packet); - void *new_data; - u32 w,h, nbelements; - int new_p; - int encoding = image->encoding; - u32 i; - - if (!pdp_packet_image_isvalid(packet)) return -1; - w = image->width; - h = image->height; - nbelements = w*h; - - new_p = pdp_packet_new_matrix(h,w, type); - if (-1 == new_p) return -1; - new_data = pdp_packet_data(new_p); - - switch (encoding){ - case PDP_IMAGE_YV12: - case PDP_IMAGE_GREY: - case PDP_IMAGE_MCHP: - /* convert first channel */ - switch(type){ - case PDP_MATRIX_TYPE_RFLOAT: - for (i=0; i<nbelements; i++) - ((float *)new_data)[i] = ((float)data[i]) * (1.0f / (float)0x8000); break; - case PDP_MATRIX_TYPE_RDOUBLE: - for (i=0; i<nbelements; i++) - ((double *)new_data)[i] = ((double)data[i]) * (1.0f / (double)0x8000); break; - case PDP_MATRIX_TYPE_CFLOAT: - for (i=0; i<nbelements; i++){ - ((float *)new_data)[i*2] = ((float)data[i]) * (1.0f / (float)0x8000); - ((float *)new_data)[i*2+1] = 0.0f; - } - break; - case PDP_MATRIX_TYPE_CDOUBLE: - for (i=0; i<nbelements; i++){ - ((double *)new_data)[i*2] = ((double)data[i]) * (1.0f / (double)0x8000); - ((double *)new_data)[i*2+1] = 0.0; - } - break; - default: - post("_pdp_packet_convert_image_to_matrix: INTERNAL ERROR"); - } - break; - default: - pdp_packet_mark_unused(new_p); - new_p = -1; - break; - } - - return new_p; - -} - -static int _pdp_packet_convert_image_to_floatmatrix(int packet, t_pdp_symbol *dest_template){ - return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RFLOAT);} -static int _pdp_packet_convert_image_to_doublematrix(int packet, t_pdp_symbol *dest_template){ - return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RDOUBLE);} -static int _pdp_packet_convert_image_to_complexfloatmatrix(int packet, t_pdp_symbol *dest_template){ - return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CFLOAT);} -static int _pdp_packet_convert_image_to_complexdoublematrix(int packet, t_pdp_symbol *dest_template){ - return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CDOUBLE);} - - -static int _pdp_packet_matrix_convert_fallback(int packet, t_pdp_symbol *dest_template) -{ - post("can't convert image type %s to %s", - pdp_packet_get_description(packet)->s_name, dest_template->s_name); - - return -1; -} - - -/* this seems like a nice spot to place a dummy gsl signal handler */ -static void _gsl_error_handler (const char * reason, - const char * file, - int line, - int gsl_errno) -{ - //post("gsl error:\nREASON: %s\nFILE:%s\nLINE:%d\nERRNO:%d", reason, file, line, gsl_errno); -} - -static int pdp_matrix_factory(t_pdp_symbol *type) -{ - return -1; -} -void pdp_matrix_setup(void) -{ - t_pdp_conversion_program *program; - - - /* setup the class */ - matrix_class = pdp_class_new(pdp_gensym("matrix/*/*/*"), pdp_matrix_factory); - matrix_class->copy = _pdp_matrix_copy; - matrix_class->clone = _pdp_matrix_clone; - matrix_class->reinit = _pdp_matrix_reinit; - matrix_class->cleanup = _pdp_matrix_cleanup; - - - /* image -> matrix: default = double/real (most ops available) */ - program = pdp_conversion_program_new(_pdp_packet_convert_image_to_doublematrix, 0); - pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/real/*"), program); - program = pdp_conversion_program_new(_pdp_packet_convert_image_to_floatmatrix, 0); - pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/real/*"), program); - program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexdoublematrix, 0); - pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/complex/*"), program); - program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexfloatmatrix, 0); - pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/complex/*"), program); - - /* matrix -> image */ - program = pdp_conversion_program_new(_pdp_packet_convert_matrix_to_greyimage, 0); - pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/grey/*"), program); - - /* fallbacks */ - program = pdp_conversion_program_new(_pdp_packet_matrix_convert_fallback, 0); - pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/*/*"), program); - pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("bitmap/*/*"), program); - pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/*/*/*/*"), program); - pdp_type_register_conversion(pdp_gensym("bitmap/*/*"), pdp_gensym("matrix/*/*/*/*"), program); - - /* setup gsl handler */ - if(gsl_set_error_handler(_gsl_error_handler)){ - post("pdp_matrix_setup: WARNING: overriding gsl error handler."); - } - -} |