aboutsummaryrefslogtreecommitdiff
path: root/system/type/pdp_matrix.c
diff options
context:
space:
mode:
Diffstat (limited to 'system/type/pdp_matrix.c')
-rw-r--r--system/type/pdp_matrix.c654
1 files changed, 654 insertions, 0 deletions
diff --git a/system/type/pdp_matrix.c b/system/type/pdp_matrix.c
new file mode 100644
index 0000000..a6fc827
--- /dev/null
+++ b/system/type/pdp_matrix.c
@@ -0,0 +1,654 @@
+
+#include <string.h>
+#include "pdp_matrix.h"
+#include "pdp_packet.h"
+#include "pdp_symbol.h"
+#include "pdp_internals.h" // for pdp_packet_new, which is not a proper constructor
+#include "pdp_post.h"
+#include "pdp_type.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:
+ pdp_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)
+{
+ pdp_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)
+{
+ //pdp_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; // is now solved through (symbol based) constructor
+ matrix_class->wakeup = _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)){
+ pdp_post("pdp_matrix_setup: WARNING: overriding gsl error handler.");
+ }
+
+}