From 213978fa8a868661dda88e6b4e7cacec1c90677a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pablo=20Mart=C3=ADn?= Date: Sun, 7 Sep 2003 21:45:05 +0000 Subject: Updating pdp to current version 0.12.2 svn path=/trunk/externals/pdp/; revision=941 --- modules/matrix_basic/Makefile | 12 + modules/matrix_basic/README | 5 + modules/matrix_basic/clusterstuff.c | 540 ++++++++++++++++++++++++++++++++++++ modules/matrix_basic/pdp_mat_lu.c | 143 ++++++++++ modules/matrix_basic/pdp_mat_mul.c | 308 ++++++++++++++++++++ modules/matrix_basic/pdp_mat_vec.c | 213 ++++++++++++++ 6 files changed, 1221 insertions(+) create mode 100644 modules/matrix_basic/Makefile create mode 100644 modules/matrix_basic/README create mode 100644 modules/matrix_basic/clusterstuff.c create mode 100644 modules/matrix_basic/pdp_mat_lu.c create mode 100644 modules/matrix_basic/pdp_mat_mul.c create mode 100644 modules/matrix_basic/pdp_mat_vec.c (limited to 'modules/matrix_basic') diff --git a/modules/matrix_basic/Makefile b/modules/matrix_basic/Makefile new file mode 100644 index 0000000..3604041 --- /dev/null +++ b/modules/matrix_basic/Makefile @@ -0,0 +1,12 @@ +current: all_modules + +include ../../Makefile.config + +PDP_MOD = pdp_mat_mul.o pdp_mat_lu.o pdp_mat_vec.o + +all_modules: $(PDP_MOD) + +clean: + rm -f *~ + rm -f *.o + diff --git a/modules/matrix_basic/README b/modules/matrix_basic/README new file mode 100644 index 0000000..d6ece1b --- /dev/null +++ b/modules/matrix_basic/README @@ -0,0 +1,5 @@ +This directory contains "normal" matrix packet processors, +derived from the t_pdp_base class defined in pdp_base.h + +Most modules are wrappers around Gnu Scientific Library (gsl) calls + diff --git a/modules/matrix_basic/clusterstuff.c b/modules/matrix_basic/clusterstuff.c new file mode 100644 index 0000000..432a807 --- /dev/null +++ b/modules/matrix_basic/clusterstuff.c @@ -0,0 +1,540 @@ +/* + * Pure Data Packet module. + * Copyright (c) 2003 by Tom Schouten + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + */ + + + +#include "pdp.h" +#include "pdp_base.h" +#include + +struct _pdp_histo; +typedef void (*t_histo_proc)(struct _pdp_histo *); + + +/* the cluster struct */ +typedef struct _cluster +{ + float N; + float cx; + float cy; +} t_cluster; + + + +typedef struct _pdp_histo +{ + t_object x_obj; + t_int x_logN; + t_symbol *x_array_sym; + t_float x_scale; + t_int x_debug; + t_int x_sample_size; /* pointcloud size */ + t_int x_nb_clusters; /* nb of clusters */ + t_cluster *x_cluster; /* cluster data (for tracking) */ + t_histo_proc x_process_method; /* what to do with the histogram */ + t_outlet *x_outlet0; + t_outlet *x_outlet1; + + /* the packet */ + int x_packet0; + + /* packet data */ + short int *x_data; + int x_width; + int x_height; + int x_nb_pixels; + + /* histo data for processor: these are stored on the stack */ + int *x_histo; + int *x_pixel_offset; + +} t_pdp_histo; + + +// join 2 clusters. clear the second one +static void cluster_join(t_cluster *cA, t_cluster *cB) +{ + float scale = 1.0f / (cA->N + cB->N); + cA->cx = (cA->N * cA->cx + cB->N * cB->cx) * scale; + cA->cy = (cA->N * cA->cy + cB->N * cB->cy) * scale; + cA->N += cB->N; + + cB->N = 0.0f; +} + +static void cluster_copy(t_cluster *cA, t_cluster *cB) +{ + cA->cx = cB->cx; + cA->cy = cB->cy; + cA->N = cB->N; +} + +static void cluster_clear(t_cluster *c) +{ + c->N = 0.0f; +} + +static void cluster_new(t_cluster *c, float x, float y) +{ + c->N = 1.0f; + c->cx = x; + c->cy = y; +} + +static float cluster_dsquared(t_cluster *cA, t_cluster *cB) +{ + float dx = cA->cx - cB->cx; + float dy = cA->cy - cB->cy; + return dx*dx + dy*dy; +} + + +static int round_up_2log(int i) +{ + int l = 0; + i--; + while (i) { + i >>= 1; + l++; + } + //post("log is %d, 2^n is %d", l, 1 << l); + + l = (l < 16) ? l : 15; + return l; +} + + +static void compute_clusters(t_pdp_histo *x) +{ + t_cluster c[x->x_sample_size]; + int i; + float scalex = 1.0f / (float)(x->x_width); + float scaley = 1.0f / (float)(x->x_height); + int nb_clusters = x->x_sample_size; + + /* build the cluster data struct */ + for (i=0; ix_sample_size; i++) + cluster_new(c+i, + ((float)(x->x_pixel_offset[i] % x->x_width)) * scalex, + ((float)(x->x_pixel_offset[i] / x->x_width)) * scaley); + + /* the clustering loop */ + while (nb_clusters > x->x_nb_clusters){ + /* initialize cA, cB, d */ + int cA=0; + int cB=1; + float d = cluster_dsquared(c+0, c+1); + int i,j; + + /* find the closest 2 clusters: + scan the distance matrix above the diagonal */ + for (i=2; ix_cluster){ + int size = sizeof(t_cluster) * x->x_nb_clusters; + x->x_cluster = (t_cluster *)pdp_alloc(size); + memcpy(x->x_cluster, c, size); + } + /* or perform tracking */ + else{ + int i,j; + /* find best matches for the first couple of clusters */ + for (i=0; ix_nb_clusters - 1; i++){ + int closest = 0; + float d_min = cluster_dsquared(x->x_cluster+i, c); + + /* get closest cluster */ + for (j=1; jx_cluster+i, c+j); + if (dj < d_min){ + closest = j; + d_min = dj; + } + } + + /* replace reference cluster with closest match */ + cluster_copy(x->x_cluster+i, c+closest); + + /* shrink matrix (like above) */ + nb_clusters--; + cluster_copy(c+closest, c+nb_clusters); + + } + /* copy the last cluster */ + cluster_copy(x->x_cluster + x->x_nb_clusters - 1, c); + } + + /* print the clusters */ + post("clusters:"); + post("\tN\tcx\tcy"); + for (i=0; ix_nb_clusters; i++){ + post("\t%d\t%0.2f\t%0.2f", + (int)x->x_cluster[i].N, + x->x_cluster[i].cx, + x->x_cluster[i].cy); + } + + + +} + +static void dump_to_array(t_pdp_histo *x) +{ + float *vec; + int nbpoints; + t_garray *a; + int i; + int *histo = x->x_histo; + int N = 1 << (x->x_logN); + float scale = 1.0f / (float)(x->x_nb_pixels); + + + /* dump to array if possible */ + if (!x->x_array_sym){ + } + + /* check if array is valid */ + else if (!(a = (t_garray *)pd_findbyclass(x->x_array_sym, garray_class))){ + post("pdp_histo: %s: no such array", x->x_array_sym->s_name); + } + /* get data */ + else if (!garray_getfloatarray(a, &nbpoints, &vec)){ + post("pdp_histo: %s: bad template", x->x_array_sym->s_name); + } + /* scale and dump in array */ + else{ + + N = (nbpoints < N) ? nbpoints : N; + for (i=0; ix_scale; + //garray_redraw(a); + } + +} + +static void get_sampleset(t_pdp_histo *x, int log_tmp_size, int threshold) +{ + int N = 1 << log_tmp_size; + int mask = N-1; + int index, nbpoints, i; + t_atom a[2]; + float scalex = 1.0f / (float)(x->x_width); + float scaley = 1.0f / (float)(x->x_height); + t_symbol *s = gensym("list"); + + /* store the offsets of the points in a in an oversized array + the oversizing is to eliminate a division and to limit the + searching for a free location after a random index is generated */ + + int offset[N]; + + /* float versions of the coordinates */ + float fx[x->x_sample_size]; + float fy[x->x_sample_size]; + float max_x, min_x, max_y, min_y; + + /* reset the array */ + memset(offset, -1, N * sizeof(int)); + + /* get the coordinates of the tempsize brightest points + and store them in a random location in the hash */ + for (i=0; ix_nb_pixels; i++){ + if (x->x_data[i] >= threshold){ + /* get a random index */ + int ri = random(); + //int ri = 0; + /* find an empty spot to store it */ + while (-1 != offset[ri & mask]) ri++; + offset[ri & mask] = i; + } + } + + + /* repack the array to get the requested + sample size at the start */ + index = 0; + nbpoints = 0; + while (nbpoints < x->x_sample_size){ + while (-1 == offset[index]) index++; // ffwd to next nonepty slot + offset[nbpoints++] = offset[index++]; // move slot + } + + /* mark output packet samples */ + memset(x->x_data, 0, 2*x->x_nb_pixels); + for (i=0; ix_sample_size; i++){ + x->x_data[offset[i]] = 0x7fff; + } + + /* send packet to left outlet */ + pdp_pass_if_valid(x->x_outlet0, &x->x_packet0); + + + /* run the clustering algo */ + x->x_pixel_offset = offset; + compute_clusters(x); + + +} + +static void get_brightest(t_pdp_histo *x) +{ + int i; + int *histo = x->x_histo; + int N = 1 << (x->x_logN); + + int index, nsamps; + + /* check requested size */ + if (x->x_sample_size > x->x_nb_pixels){ + post("WARNING: more samples requested than pixels in image"); + x->x_sample_size = x->x_nb_pixels; + } + + + /* find limiting index */ + index = N; + nsamps = 0; + while (nsamps < x->x_sample_size){ + index--; + nsamps += histo[index]; + } + + /* status report */ + if (x->x_debug){ + post("found %d samples between h[%d] and h[%d]", nsamps, index, N-1); + } + + /* get a representative set from the candidates + the tempbuf is the rounded log of the nb of samples + 1 + so it is at least 50% sparse */ + get_sampleset(x, round_up_2log(nsamps) + 1, index << (15-x->x_logN)); + +} + + +static void _pdp_histo_perform(t_pdp_histo *x) +{ + short int *pp; + int N = 1 << x->x_logN; + int nbpixels = x->x_width * x->x_height, i; + + int histo[N]; + + /* init */ + for (i=0; ix_data[i] >> (15 - x->x_logN); + if (index < 0) index = 0; /* negative -> zero */ + histo[index]++; + } + + /* save the histo stack location */ + x->x_histo = histo; + + /* print it */ + if (x->x_debug){ + post("histogram:"); + for (i=0; ix_process_method(x); + + +} + + +// packet is an image/*/* packet or invalid */ +static void pdp_histo_perform(t_pdp_histo *x) +{ + t_pdp *header0 = pdp_packet_header(x->x_packet0); + void *data0 = pdp_packet_data(x->x_packet0); + if (!header0 || !data0) return; + + x->x_width = header0->info.image.width; + x->x_height = header0->info.image.height; + x->x_nb_pixels = x->x_width * x->x_height; + x->x_data = data0; + + _pdp_histo_perform(x); +} + + + +static void pdp_histo_input_0(t_pdp_histo *x, t_symbol *s, t_floatarg f) +{ + int packet = (int)f; + + /* register */ + if (s == gensym("register_ro")){ + /* replace if not compatible or we are not interpolating */ + pdp_packet_mark_unused(x->x_packet0); + x->x_packet0 = pdp_packet_convert_rw(packet, pdp_gensym("image/grey/*")); + + } + + if (s == gensym("process")){ + pdp_histo_perform(x); + } + +} + + + +static void pdp_histo_samplesize(t_pdp_histo *x, t_floatarg f) +{ + int i = (int)f; + if (i >= x->x_nb_clusters ) x->x_sample_size = i; +} + +static void pdp_histo_clusters(t_pdp_histo *x, t_floatarg f) +{ + int i = (int)f; + if (i>=2 && i<= x->x_sample_size){ + x->x_nb_clusters = i; + if (x->x_cluster) pdp_dealloc(x->x_cluster); + x->x_cluster = 0; + } +} +static void pdp_histo_scale(t_pdp_histo *x, t_floatarg f){x->x_scale = f;} + + + +static void pdp_histo_size(t_pdp_histo *x, t_floatarg f) +{ + int i = (int)f; + if (i < 1) return; + x->x_logN = round_up_2log(i); +} + + +static void pdp_histo_array(t_pdp_histo *x, t_symbol *s) +{ + //post("setting symbol %x", s); + x->x_array_sym = s; +} + + +static void pdp_histo_free(t_pdp_histo *x) +{ + pdp_packet_mark_unused(x->x_packet0); + if (x->x_cluster) pdp_dealloc(x->x_cluster); +} + + +t_class *pdp_histo_class; + + + +void *pdp_histo_new(t_floatarg f) +{ + + t_pdp_histo *x = (t_pdp_histo *)pd_new(pdp_histo_class); + if (f == 0.0f) f = 64; + pdp_histo_size(x, f); + x->x_packet0 = -1; + x->x_debug = 0; + x->x_sample_size = 16; + x->x_nb_clusters = 3; + x->x_cluster = 0; + return (void *)x; +} + + +void *pdp_histo_array_new(t_symbol *s, t_float f, t_float f2) +{ + t_pdp_histo *x = (t_pdp_histo *)pdp_histo_new(f); + if (f2 == 0.0f) f2 = 1.0f; + pdp_histo_scale(x, f2); + pdp_histo_array(x, s); + x->x_process_method = dump_to_array; + return (void *)x; +} + +void *pdp_histo_sample_new(t_float nbsamples, t_float histosize) +{ + t_pdp_histo *x; + if (histosize == 0.0f) histosize = 256.0f; + x = (t_pdp_histo *)pdp_histo_new(histosize); + if (nbsamples == 0.0f) nbsamples = 16.0f; + pdp_histo_samplesize(x, nbsamples); + x->x_process_method = get_brightest; + x->x_outlet0 = outlet_new(&x->x_obj, gensym("anything")); + //x->x_outlet1 = outlet_new(&x->x_obj, gensym("anything")); + + inlet_new((t_object *)x, (t_pd *)&x->x_obj, gensym("float"), gensym("nbpoints")); + + return (void *)x; +} + + +#ifdef __cplusplus +extern "C" +{ +#endif + + +void pdp_histo_setup(void) +{ + + pdp_histo_class = class_new(gensym("pdp_histo"), (t_newmethod)pdp_histo_array_new, + (t_method)pdp_histo_free, sizeof(t_pdp_histo), 0, A_DEFSYMBOL, A_DEFFLOAT, A_DEFFLOAT, A_NULL); + + class_addcreator((t_newmethod)pdp_histo_sample_new, gensym("pdp_pointcloud"), A_DEFFLOAT, A_DEFFLOAT, A_NULL); + + class_addmethod(pdp_histo_class, (t_method)pdp_histo_input_0, gensym("pdp"), A_SYMBOL, A_DEFFLOAT, A_NULL); + + class_addmethod(pdp_histo_class, (t_method)pdp_histo_size, gensym("size"), A_FLOAT, A_NULL); + class_addmethod(pdp_histo_class, (t_method)pdp_histo_size, gensym("scale"), A_FLOAT, A_NULL); + class_addmethod(pdp_histo_class, (t_method)pdp_histo_array, gensym("array"), A_SYMBOL, A_NULL); + class_addmethod(pdp_histo_class, (t_method)pdp_histo_samplesize, gensym("nbpoints"), A_FLOAT, A_NULL); + class_addmethod(pdp_histo_class, (t_method)pdp_histo_clusters, gensym("nbclusters"), A_FLOAT, A_NULL); +} + +#ifdef __cplusplus +} +#endif + + + diff --git a/modules/matrix_basic/pdp_mat_lu.c b/modules/matrix_basic/pdp_mat_lu.c new file mode 100644 index 0000000..af8931d --- /dev/null +++ b/modules/matrix_basic/pdp_mat_lu.c @@ -0,0 +1,143 @@ +/* + * Pure Data Packet module. LU decomposition module. + * Copyright (c) by Tom Schouten + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + */ + +//#include +//#include +//#include +//#include +#include "pdp.h" +#include "pdp_base.h" + + +typedef struct pdp_mat_LU_struct +{ + t_pdp_base x_base; + +} t_pdp_mat_LU; + + + +static void pdp_mat_LU_process_LU_inverse(t_pdp_mat_LU *x) +{ + int p = pdp_base_get_packet(x, 0); + int p_LU = pdp_packet_matrix_LU_to_inverse(p); + pdp_base_set_packet(x, 0, p_LU); // replace packet +} + +static void pdp_mat_LU_process_LU(t_pdp_mat_LU *x) +{ + int p = pdp_base_get_packet(x, 0); + int p_LU = pdp_packet_matrix_LU(p); + pdp_base_set_packet(x, 0, p_LU); // replace packet +} + +static void pdp_mat_LU_process_LU_solve(t_pdp_mat_LU *x) +{ + int p0 = pdp_base_get_packet(x, 0); + int p1 = pdp_base_get_packet(x, 1); + int pvr, pm, pv; + + /* determine which is vector and which is matrix */ + if (pdp_packet_matrix_ismatrix(p0) && pdp_packet_matrix_isvector(p1)){ + pm = p0; + pv = p1; + } + else { + pm = p1; + pv = p0; + } + + /* create the result vector */ + pvr = pdp_packet_matrix_LU_solve(pm, pv); + + /* replace the active packet */ + pdp_base_set_packet(x, 0, pvr); +} + +static void pdp_mat_LU_free(t_pdp_mat_LU *x) +{ + /* remove process method from queue before deleting data */ + pdp_base_free(x); +} + +t_class *pdp_mat_LU_class; + + +/* common new methods */ +t_pdp_mat_LU *pdp_mat_LU_base_new(void) +{ + t_pdp_mat_LU *x = (t_pdp_mat_LU *)pd_new(pdp_mat_LU_class); + pdp_base_init(x); + pdp_base_add_pdp_outlet(x); + return x; +} + +void *pdp_mat_LU_inverse_new(void) +{ + t_pdp_mat_LU *x = pdp_mat_LU_base_new(); + pdp_base_set_process_method(x,(t_pdp_method)pdp_mat_LU_process_LU_inverse); + pdp_base_readonly_active_inlet(x); + return (void *)x; +} + + +void *pdp_mat_LU_new(void) +{ + t_pdp_mat_LU *x = pdp_mat_LU_base_new(); + pdp_base_set_process_method(x,(t_pdp_method)pdp_mat_LU_process_LU); + pdp_base_readonly_active_inlet(x); + return (void *)x; +} + +void *pdp_mat_LU_solve_new(void) +{ + t_pdp_mat_LU *x = pdp_mat_LU_base_new(); + pdp_base_set_process_method(x,(t_pdp_method)pdp_mat_LU_process_LU_solve); + pdp_base_readonly_active_inlet(x); + pdp_base_add_pdp_inlet(x); + return (void *)x; +} + +#ifdef __cplusplus +extern "C" +{ +#endif + + +void pdp_mat_lu_setup(void) +{ + + + pdp_mat_LU_class = class_new(gensym("pdp_m_LU_inverse"), (t_newmethod)pdp_mat_LU_inverse_new, + (t_method)pdp_mat_LU_free, sizeof(t_pdp_mat_LU), 0, A_NULL); + + pdp_base_setup(pdp_mat_LU_class); + + + class_addcreator((t_newmethod)pdp_mat_LU_new, gensym("pdp_m_LU"), A_NULL); + class_addcreator((t_newmethod)pdp_mat_LU_solve_new, gensym("pdp_m_LU_solve"), A_NULL); + + + +} + +#ifdef __cplusplus +} +#endif diff --git a/modules/matrix_basic/pdp_mat_mul.c b/modules/matrix_basic/pdp_mat_mul.c new file mode 100644 index 0000000..329813f --- /dev/null +++ b/modules/matrix_basic/pdp_mat_mul.c @@ -0,0 +1,308 @@ +/* + * Pure Data Packet module. Matrix multiplication module + * Copyright (c) by Tom Schouten + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + */ + +//#include +//#include +//#include +//#include +#include "pdp.h" +#include "pdp_base.h" + + +typedef struct pdp_mat_mm_struct +{ + t_pdp_base x_base; + CBLAS_TRANSPOSE_t x_T0; + CBLAS_TRANSPOSE_t x_T1; + int x_M0; + int x_M1; + + float x_scale_r; + float x_scale_i; + +} t_pdp_mat_mm; + + +static void pdp_mat_mm_rscale(t_pdp_mat_mm *x, t_floatarg r) +{ + x->x_scale_r = r; + x->x_scale_i = 0.0f; +} + +static void pdp_mat_mm_cscale(t_pdp_mat_mm *x, t_floatarg r, t_floatarg i) +{ + x->x_scale_r = r; + x->x_scale_i = i; +} + + +/* matrix multilpy */ +static void pdp_mat_mv_process_mul(t_pdp_mat_mm *x) +{ + int pA = pdp_base_get_packet(x, 0); + int pB = pdp_base_get_packet(x, 1); + int p0, p1, pR; + + /* determine which one is the vector */ + if (pdp_packet_matrix_isvector(pA)){ + p0 = pB; + p1 = pA; + } + else { + p1 = pB; + p0 = pA; + } + + pR = pdp_packet_new_matrix_product_result(x->x_T0, CblasNoTrans, p0, p1); + + if (-1 != pR){ + pdp_packet_matrix_setzero(pR); + if (pdp_packet_matrix_blas_mv(x->x_T0, p0, p1, pR, x->x_scale_r, x->x_scale_i)){ + //post("pdp_packet_matrix_blas_mm failed"); + pdp_packet_mark_unused(pR); + pR = -1; + } + } + else { + //post("pdp_packet_new_matrix_product_result failed"); + } + + /* replace with result */ + pdp_base_set_packet(x, 0, pR); + +} + +/* matrix vector multilpy */ +static void pdp_mat_mm_process_mul(t_pdp_mat_mm *x) +{ + int pA = pdp_base_get_packet(x, 0); + int pB = pdp_base_get_packet(x, 1); + int p0, p1, pR; + + p0 = (x->x_M0) ? pB : pA; + p1 = (x->x_M1) ? pB : pA; + + pR = pdp_packet_new_matrix_product_result(x->x_T0, x->x_T1, p0, p1); + + if (-1 != pR){ + pdp_packet_matrix_setzero(pR); + if (pdp_packet_matrix_blas_mm(x->x_T0, x->x_T1, p0, p1, pR, x->x_scale_r, x->x_scale_i)){ + //post("pdp_packet_matrix_blas_mm failed"); + pdp_packet_mark_unused(pR); + pR = -1; + } + } + else { + //post("pdp_packet_new_matrix_product_result failed"); + } + + /* replace with result */ + pdp_base_set_packet(x, 0, pR); + +} +/* matrix macc */ +static void pdp_mat_mm_process_mac(t_pdp_mat_mm *x) +{ + int pC = pdp_base_get_packet(x, 0); + int pA = pdp_base_get_packet(x, 1); + int pB = pdp_base_get_packet(x, 2); + int p0, p1; + + p0 = (x->x_M0) ? pB : pA; + p1 = (x->x_M1) ? pB : pA; + + if (pdp_packet_matrix_blas_mm(x->x_T0, x->x_T1, p0, p1, pC, x->x_scale_r, x->x_scale_i)){ + //post("pdp_packet_matrix_blas_mm failed"); + pdp_base_set_packet(x, 0, -1); // delete packet + } + +} + + +static void pdp_mat_mm_free(t_pdp_mat_mm *x) +{ + /* remove process method from queue before deleting data */ + pdp_base_free(x); +} + +t_class *pdp_mat_mm_class; + + +/* common new method */ +void *pdp_mat_mm_new(void) +{ + int i; + t_pdp_mat_mm *x = (t_pdp_mat_mm *)pd_new(pdp_mat_mm_class); + + /* super init */ + pdp_base_init(x); + + /* outlet */ + pdp_base_add_pdp_outlet(x); + + + return (void *)x; +} + + +static int pdp_mat_mm_setup_routing_M0(t_pdp_mat_mm *x, t_symbol *s0) +{ + if ('A' == s0->s_name[0]){x->x_M0 = 0;} else if ('B' == s0->s_name[0]) {x->x_M0 = 1;} else return 0; + + if ((gensym("A") == s0) || (gensym("B") == s0)) x->x_T0 = CblasNoTrans; + else if ((gensym("A^T") == s0) || (gensym("B^T") == s0)) x->x_T0 = CblasConjTrans; + else if ((gensym("A^H") == s0) || (gensym("B^H") == s0)) x->x_T0 = CblasConjTrans; + else return 0; + + return 1; +} + +static int pdp_mat_mm_setup_routing_M1(t_pdp_mat_mm *x, t_symbol *s1) +{ + + if ('A' == s1->s_name[0]){x->x_M1 = 0;} else if ('B' == s1->s_name[0]) {x->x_M1 = 1;} else return 0; + + /* setup second matrix transpose operation */ + if ((gensym("A") == s1) || (gensym("B") == s1)) x->x_T1 = CblasNoTrans; + else if ((gensym("A^T") == s1) || (gensym("B^T") == s1)) x->x_T1 = CblasConjTrans; + else if ((gensym("A^H") == s1) || (gensym("B^H") == s1)) x->x_T1 = CblasConjTrans; + else return 0; + + return 1; +} + + +static int pdp_mat_mm_setup_scaling(t_pdp_mat_mm *x, t_symbol *scale) +{ + int success = 1; + + /* setup scaling inlet */ + if ((gensym ("rscale") == scale) || (gensym("r") == scale)){ + pdp_base_add_gen_inlet(x, gensym("float"), gensym("rscale")); + } + else if ((gensym ("cscale") == scale) || (gensym("c") == scale)){ + pdp_base_add_gen_inlet(x, gensym("list"), gensym("cscale")); + } + else if (gensym ("") != scale) success = 0; + + return success; +} + +void *pdp_mat_mm_new_mul_common(t_symbol *s0, t_symbol *s1, t_symbol *scale, int ein) +{ + t_pdp_mat_mm *x = pdp_mat_mm_new(); + + /* add extra pdp inlets */ + while (ein--) pdp_base_add_pdp_inlet(x); + + /* setup routing */ + if (!pdp_mat_mm_setup_routing_M0(x, s0)) goto error; + if (!pdp_mat_mm_setup_routing_M1(x, s1)) goto error; + if (!pdp_mat_mm_setup_scaling(x, scale)) goto error; + + /* default scale = 1 */ + pdp_mat_mm_cscale(x, 1.0f, 0.0f); + return (void *)x; + + error: + pd_free((void *)x); + return 0; +} + +void *pdp_mat_mv_new_mul_common(t_symbol *s0, t_symbol *scale, int ein) +{ + t_pdp_mat_mm *x = pdp_mat_mm_new(); + + /* add extra pdp inlets */ + while (ein--) pdp_base_add_pdp_inlet(x); + + /* setup routing */ + if (!pdp_mat_mm_setup_routing_M0(x, s0)) goto error; + if (!pdp_mat_mm_setup_scaling(x, scale)) goto error; + + /* default scale = 1 */ + pdp_mat_mm_cscale(x, 1.0f, 0.0f); + return (void *)x; + + error: + pd_free((void *)x); + return 0; +} + +void *pdp_mat_mm_new_mul(t_symbol *s0, t_symbol *s1, t_symbol *scale) +{ + t_pdp_mat_mm *x = pdp_mat_mm_new_mul_common(s0, s1, scale, 1); + if(x){ + pdp_base_set_process_method(x, (t_pdp_method)pdp_mat_mm_process_mul); + pdp_base_readonly_active_inlet(x); + } + return x; +} + +void *pdp_mat_mv_new_mul(t_symbol *s0, t_symbol *scale) +{ + t_pdp_mat_mm *x = pdp_mat_mv_new_mul_common(s0, scale, 1); + if(x){ + pdp_base_set_process_method(x, (t_pdp_method)pdp_mat_mv_process_mul); + pdp_base_readonly_active_inlet(x); + } + return x; +} + +void *pdp_mat_mm_new_mac(t_symbol *s0, t_symbol *s1, t_symbol *scale) +{ + t_pdp_mat_mm *x = pdp_mat_mm_new_mul_common(s0, s1, scale, 2); + if (x){ + pdp_base_set_process_method(x, (t_pdp_method)pdp_mat_mm_process_mac); + } + return x; +} + + +#ifdef __cplusplus +extern "C" +{ +#endif + + +void pdp_mat_mul_setup(void) +{ + + + pdp_mat_mm_class = class_new(gensym("pdp_m_mm"), (t_newmethod)pdp_mat_mm_new_mul, + (t_method)pdp_mat_mm_free, sizeof(t_pdp_mat_mm), 0, A_SYMBOL, A_SYMBOL, A_DEFSYMBOL, A_NULL); + + pdp_base_setup(pdp_mat_mm_class); + + class_addcreator((t_newmethod)pdp_mat_mm_new_mac, gensym("pdp_m_+=mm"), + A_SYMBOL, A_SYMBOL, A_DEFSYMBOL, A_NULL); + + class_addcreator((t_newmethod)pdp_mat_mv_new_mul, gensym("pdp_m_mv"), + A_SYMBOL, A_DEFSYMBOL, A_NULL); + + + class_addmethod(pdp_mat_mm_class, (t_method)pdp_mat_mm_rscale, gensym("rscale"), A_FLOAT, A_NULL); + class_addmethod(pdp_mat_mm_class, (t_method)pdp_mat_mm_cscale, gensym("cscale"), A_FLOAT, A_FLOAT, A_NULL); + +} + +#ifdef __cplusplus +} +#endif diff --git a/modules/matrix_basic/pdp_mat_vec.c b/modules/matrix_basic/pdp_mat_vec.c new file mode 100644 index 0000000..16257c8 --- /dev/null +++ b/modules/matrix_basic/pdp_mat_vec.c @@ -0,0 +1,213 @@ +/* + * Pure Data Packet module. Vector modules. + * Copyright (c) by Tom Schouten + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + */ + +//#include +//#include +//#include +//#include +#include "pdp.h" +#include "pdp_base.h" + + +typedef struct pdp_mat_vec_struct +{ + t_pdp_base x_base; + int x_type; + t_outlet *x_out; + int x_accept_list; + int x_list_size; + t_atom *x_list; + +} t_pdp_mat_vec; + + +#define GETFLOAT(x) ((x)->a_type == A_FLOAT ? (x)->a_w.w_float : 0.0f) +#define GETDOUBLE(x) (double)GETFLOAT(x) + + +static void pdp_mat_vec_list_in(t_pdp_mat_vec *x, t_symbol *s, int argc, t_atom *argv) +{ + int i; + int vp = -1; + int f; + int dim = argc; + + if (!x->x_accept_list) return; //check if this handler is enabled + if (!argc) return; //reject empty list + + switch(x->x_type){ + case PDP_MATRIX_TYPE_CFLOAT: + if (argc & 1) return; //reject odd nb elements + dim >>= 1; //halve dimension + case PDP_MATRIX_TYPE_RFLOAT: + vp = pdp_packet_new_matrix(dim, 1, x->x_type); + if (-1 != vp){ + float *data = (float *)pdp_packet_data(vp); + for (i=0; i>= 1; //halve dimension + case PDP_MATRIX_TYPE_RDOUBLE: + vp = pdp_packet_new_matrix(dim, 1, x->x_type); + if (-1 != vp){ + double *data = (double *)pdp_packet_data(vp); + for (i=0; isize; + if ((type == PDP_MATRIX_TYPE_CFLOAT) + || (type == PDP_MATRIX_TYPE_CDOUBLE)) + outlist_size <<= 1; + + /* realloc list if necessary */ + if (outlist_size > x->x_list_size){ + free(x->x_list); + x->x_list = (t_atom *)pdp_alloc(sizeof(t_atom) * outlist_size); + x->x_list_size = outlist_size; + } + + /* copy data */ + switch(type){ + case PDP_MATRIX_TYPE_RFLOAT: + case PDP_MATRIX_TYPE_CFLOAT: + fdata = (float *)pdp_packet_data(p); + for (i=0; ix_list[i], fdata[i]); + break; + + case PDP_MATRIX_TYPE_RDOUBLE: + case PDP_MATRIX_TYPE_CDOUBLE: + ddata = (double *)pdp_packet_data(p); + for (i=0; ix_list[i], (float)ddata[i]); + break; + + } + + /* dispose of vector packet and output list */ + pdp_packet_mark_unused(p); + outlet_list(x->x_out, &s_list, outlist_size, x->x_list); +} + +static void pdp_mat_vec_free(t_pdp_mat_vec *x) +{ + /* remove process method from queue before deleting data */ + pdp_base_free(x); + + /* delete list */ + if (x->x_list) pdp_dealloc (x->x_list); +} + +t_class *pdp_mat_vec_class; + + +/* common new methods */ +t_pdp_mat_vec *pdp_mat_vec_base_new(void) +{ + t_pdp_mat_vec *x = (t_pdp_mat_vec *)pd_new(pdp_mat_vec_class); + pdp_base_init(x); + x->x_type = PDP_MATRIX_TYPE_CFLOAT; + x->x_accept_list = 0; + x->x_list_size = 0; + x->x_list = 0; + return x; +} + +void *pdp_mat_vec_list2vec_new(t_symbol *type) +{ + t_pdp_mat_vec *x = pdp_mat_vec_base_new(); + pdp_base_disable_active_inlet(x); + pdp_base_add_pdp_outlet(x); + x->x_accept_list = 1; + if ((gensym ("") == type) || (gensym ("double/real") == type)) x->x_type = PDP_MATRIX_TYPE_RDOUBLE; + else if (gensym ("double/complex") == type) x->x_type = PDP_MATRIX_TYPE_CDOUBLE; + else if (gensym ("float/real") == type) x->x_type = PDP_MATRIX_TYPE_RFLOAT; + else if (gensym ("float/complex") == type) x->x_type = PDP_MATRIX_TYPE_CFLOAT; + else { + pd_free((t_pd *)x); + x = 0; + } + return (void *)x; +} + + +void *pdp_mat_vec_vec2list_new(t_symbol *type) +{ + t_pdp_mat_vec *x = pdp_mat_vec_base_new(); + x->x_out = outlet_new((t_object *)x, &s_anything); + pdp_base_set_postproc_method(x,(t_pdp_method)pdp_mat_vec_list_out); + pdp_base_readonly_active_inlet(x); + return (void *)x; +} + + +#ifdef __cplusplus +extern "C" +{ +#endif + + +void pdp_mat_vec_setup(void) +{ + + + pdp_mat_vec_class = class_new(gensym("pdp_m_list2vec"), (t_newmethod)pdp_mat_vec_list2vec_new, + (t_method)pdp_mat_vec_free, sizeof(t_pdp_mat_vec), 0, A_DEFSYMBOL, A_NULL); + + pdp_base_setup(pdp_mat_vec_class); + + + class_addcreator((t_newmethod)pdp_mat_vec_vec2list_new, gensym("pdp_m_vec2list"), A_NULL); + + class_addlist(pdp_mat_vec_class, (t_method)pdp_mat_vec_list_in); + + +} + +#ifdef __cplusplus +} +#endif -- cgit v1.2.1