From b694c274836ac8b04d644711ac324eac2e9ab83e Mon Sep 17 00:00:00 2001 From: Hans-Christoph Steiner Date: Fri, 16 Dec 2005 01:05:40 +0000 Subject: checking in pdp 0.12.4 from http://zwizwa.fartit.com/pd/pdp/pdp-0.12.4.tar.gz svn path=/trunk/externals/pdp/; revision=4232 --- system/type/pdp_matrix.c | 654 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 654 insertions(+) create mode 100644 system/type/pdp_matrix.c (limited to 'system/type/pdp_matrix.c') 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 +#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; iperm.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; iencoding; + 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; is_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."); + } + +} -- cgit v1.2.1