Updating pdp to current version 0.12.2

svn path=/trunk/externals/pdp/; revision=941

6 files changed, 1221 insertions, 0 deletions

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 <pdp@zzz.kotnet.org>

+ *

+ *   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 <math.h>

+

+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; i<x->x_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; i<nb_clusters; i++){ 

+	    for (j=0; j<i; j++){

+		float dij = cluster_dsquared(c+i, c+j);

+		if (dij < d){

+		    cA = j;

+		    cB = i;

+		    d = dij;

+		}

+	    }

+	}

+

+	/* join the two clusters (cA < cB) */

+	cluster_join (c+cA, c+cB);

+

+	/* reduce the distance matrix by moving

+	   the last element to the empty spot cB */

+	nb_clusters--;

+	cluster_copy (c+cB, c+nb_clusters);

+    }

+

+    /* copy cluster data */

+    if (!x->x_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; i<x->x_nb_clusters - 1; i++){

+	    int closest = 0;

+	    float d_min = cluster_dsquared(x->x_cluster+i, c);

+

+	    /* get closest cluster */

+	    for (j=1; j<nb_clusters; j++){

+		float dj = cluster_dsquared(x->x_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; i<x->x_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; i<N; i++) vec[i] = (float)(histo[i]) * scale * x->x_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; i<x->x_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; i<x->x_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; i<N; i++) histo[i] = 0;

+

+    /* build histo */

+    for (i=0; i<nbpixels; i++){

+	int index = x->x_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; i<N; i++){

+	    fprintf(stderr, "%d\t", histo[i]);

+	    if (!(i % 10)) post("");

+	}

+	post("");

+    }

+

+    /* call the processor */

+    x->x_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 <pdp@zzz.kotnet.org>
+ *
+ *   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 <gsl/gsl_block.h>
+//#include <gsl/gsl_vector.h>
+//#include <gsl/gsl_matrix.h>
+//#include <gsl/gsl_blas.h>
+#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 <pdp@zzz.kotnet.org>
+ *
+ *   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 <gsl/gsl_block.h>
+//#include <gsl/gsl_vector.h>
+//#include <gsl/gsl_matrix.h>
+//#include <gsl/gsl_blas.h>
+#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 <pdp@zzz.kotnet.org>
+ *
+ *   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 <gsl/gsl_block.h>
+//#include <gsl/gsl_vector.h>
+//#include <gsl/gsl_matrix.h>
+//#include <gsl/gsl_blas.h>
+#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<argc; i++) data[i] = GETFLOAT(&argv[i]);
+	}
+	break;
+    case PDP_MATRIX_TYPE_CDOUBLE:
+	if (argc & 1) return; //reject odd nb elements
+	dim >>= 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; i<argc; i++) data[i] = GETDOUBLE(&argv[i]);
+	}
+	break;
+    default:
+	break;
+    }
+
+    if (-1 != vp){
+	/* store vector packet */
+	pdp_base_set_packet(x, 0, vp);
+	pdp_base_bang(x);
+    }
+}
+
+static void pdp_mat_vec_list_out(t_pdp_mat_vec *x)
+{
+    int p = pdp_base_move_packet(x, 0);
+    int type = pdp_packet_matrix_get_type(p);
+    int outlist_size;
+    float *fdata = 0;
+    double *ddata = 0;
+    int i;
+
+    /* check if it's a vector */
+    gsl_vector *m = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p, type);
+    if (!pdp_packet_matrix_isvector(p)) return;
+
+    /* get list size */
+    outlist_size = m->size;
+    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; i<outlist_size; i++)
+	    SETFLOAT(&x->x_list[i], fdata[i]);
+	break;
+	
+    case PDP_MATRIX_TYPE_RDOUBLE:
+    case PDP_MATRIX_TYPE_CDOUBLE:
+	ddata = (double *)pdp_packet_data(p);
+	for (i=0; i<outlist_size; i++)
+	    SETFLOAT(&x->x_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