aboutsummaryrefslogtreecommitdiff
path: root/modules/matrix_basic
diff options
context:
space:
mode:
authorHans-Christoph Steiner <eighthave@users.sourceforge.net>2005-12-15 07:26:47 +0000
committerHans-Christoph Steiner <eighthave@users.sourceforge.net>2005-12-15 07:26:47 +0000
commit37b6643df2df7d784a31ca73f7bb90dc109c2401 (patch)
treea8664e5adcfcb60cae136063d627549ecb76619b /modules/matrix_basic
parentc50ce0e0217ea07e2d450add2ab29cecea66fa96 (diff)
removing PDP source (except debian files) before import of PDP 0.12.4
svn path=/trunk/externals/pdp/; revision=4217
Diffstat (limited to 'modules/matrix_basic')
-rw-r--r--modules/matrix_basic/Makefile12
-rw-r--r--modules/matrix_basic/README5
-rw-r--r--modules/matrix_basic/clusterstuff.c540
-rw-r--r--modules/matrix_basic/pdp_mat_lu.c143
-rw-r--r--modules/matrix_basic/pdp_mat_mul.c308
-rw-r--r--modules/matrix_basic/pdp_mat_vec.c213
6 files changed, 0 insertions, 1221 deletions
diff --git a/modules/matrix_basic/Makefile b/modules/matrix_basic/Makefile
deleted file mode 100644
index 3604041..0000000
--- a/modules/matrix_basic/Makefile
+++ /dev/null
@@ -1,12 +0,0 @@
-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
deleted file mode 100644
index d6ece1b..0000000
--- a/modules/matrix_basic/README
+++ /dev/null
@@ -1,5 +0,0 @@
-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
deleted file mode 100644
index 432a807..0000000
--- a/modules/matrix_basic/clusterstuff.c
+++ /dev/null
@@ -1,540 +0,0 @@
-/*
- * 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
deleted file mode 100644
index af8931d..0000000
--- a/modules/matrix_basic/pdp_mat_lu.c
+++ /dev/null
@@ -1,143 +0,0 @@
-/*
- * 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
deleted file mode 100644
index 329813f..0000000
--- a/modules/matrix_basic/pdp_mat_mul.c
+++ /dev/null
@@ -1,308 +0,0 @@
-/*
- * 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
deleted file mode 100644
index 16257c8..0000000
--- a/modules/matrix_basic/pdp_mat_vec.c
+++ /dev/null
@@ -1,213 +0,0 @@
-/*
- * 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