3
#include "pdp_matrix.h"
4
#include "pdp_packet.h"
5
#include "pdp_symbol.h"
6
#include "pdp_internals.h" // for pdp_packet_new, which is not a proper constructor
10
/* the class object */
11
static t_pdp_class *matrix_class;
15
/* declare header and subheader variables and exit with errval if invalid */
16
#define VALID_MATRIX_HEADER(packet, header, subheader, mtype, errval) \
17
t_pdp * header = pdp_packet_header( packet ); \
18
t_matrix * subheader = (t_matrix *)pdp_packet_subheader( packet ); \
19
if (! header ) return errval; \
20
if (PDP_MATRIX != header->type) return errval; \
21
if (mtype) {if (subheader->type != mtype) return errval;}
24
int pdp_packet_matrix_isvalid(int p)
26
VALID_MATRIX_HEADER(p, h, m, 0, 0);
31
void *pdp_packet_matrix_get_gsl_matrix(int p, u32 type)
33
VALID_MATRIX_HEADER(p, h, m, type, 0);
37
void *pdp_packet_matrix_get_gsl_vector(int p, u32 type)
39
VALID_MATRIX_HEADER(p, h, m, type, 0);
43
int pdp_packet_matrix_get_type(int p)
45
VALID_MATRIX_HEADER(p, h, m, 0, 0);
49
int pdp_packet_matrix_isvector(int p)
51
VALID_MATRIX_HEADER(p, h, m, 0, 0);
52
return ((m->rows == 1) || (m->columns == 1));
55
int pdp_packet_matrix_ismatrix(int p)
57
VALID_MATRIX_HEADER(p, h, m, 0, 0);
58
return ((m->rows != 1) && (m->columns != 1));
64
/* gsl blas matrix/vector multiplication:
68
(const gsl_vector * x, const gsl_vector * y, double * result)
95
CBLAS_TRANSPOSE_t TransA,
96
CBLAS_TRANSPOSE_t TransB,
111
/* compute the matrix inverse using the LU decomposition */
112
/* it only works for double real/complex */
113
int pdp_packet_matrix_LU_to_inverse(int p)
117
int type = pdp_packet_matrix_get_type(p);
118
t_matrix *sheader = (t_matrix *)pdp_packet_subheader(p);
119
gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type);
123
if (n != m->gsl_columns) return -1;
124
new_p = pdp_packet_new_matrix(n, n, type);
125
if (-1 == new_p) return -1;
126
im = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(new_p, type);
129
case PDP_MATRIX_TYPE_RDOUBLE:
130
gsl_linalg_LU_invert (m, &sheader->perm, im); break;
131
case PDP_MATRIX_TYPE_CDOUBLE:
132
gsl_linalg_complex_LU_invert ((gsl_matrix_complex *)m, &sheader->perm, (gsl_matrix_complex *)im); break;
134
pdp_packet_mark_unused(new_p);
139
/* compute the LU decomposition of a square matrix */
140
/* it only works for double real/complex */
141
int pdp_packet_matrix_LU(int p)
145
int type = pdp_packet_matrix_get_type(p);
147
t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p);
148
gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type);
152
if (n != m->gsl_columns) return -1;
153
p_LU = pdp_packet_new_matrix(n, n, type);
154
if (-1 == p_LU) return -1;
155
sh_m_LU = (t_matrix *)pdp_packet_subheader(p_LU);
156
m_LU = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_LU, type);
157
/* copy matrix data: move this to copy method */
158
memcpy(pdp_packet_data(p_LU), pdp_packet_data(p), sh_m->block.size);
161
case PDP_MATRIX_TYPE_RDOUBLE:
162
gsl_linalg_LU_decomp (m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break;
163
case PDP_MATRIX_TYPE_CDOUBLE:
164
gsl_linalg_complex_LU_decomp ((gsl_matrix_complex *)m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break;
166
pdp_packet_mark_unused(p_LU);
172
int pdp_packet_matrix_LU_solve(int p_matrix, int p_vector)
174
int type = pdp_packet_matrix_get_type(p_matrix);
175
int p_result_vector = pdp_packet_clone_rw(p_vector);
176
gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_matrix, type);
177
gsl_vector *v = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_vector, type);
178
gsl_vector *rv = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_result_vector, type);
179
t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p_matrix);
181
if (!(m && v && rv)) goto error;
184
case PDP_MATRIX_TYPE_RDOUBLE:
185
if (gsl_linalg_LU_solve (m, &sh_m->perm, v, rv)) goto error; break;
186
case PDP_MATRIX_TYPE_CDOUBLE:
187
if(gsl_linalg_complex_LU_solve ((gsl_matrix_complex*)m, &sh_m->perm,
188
(gsl_vector_complex *)v, (gsl_vector_complex *)rv)) goto error; break;
192
return p_result_vector;
195
pdp_packet_mark_unused(p_result_vector);
200
/* matrix matrix mul: C is defining type
201
returns 0 on success */
202
int pdp_packet_matrix_blas_mm
204
CBLAS_TRANSPOSE_t TransA,
205
CBLAS_TRANSPOSE_t TransB,
213
gsl_complex_float cf_scale = {{scale_r, scale_i}};
214
gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}};
215
gsl_complex_float cf_one = {{1.0f, 0.0f}};
216
gsl_complex cd_one = {{1.0, 0.0}};
217
gsl_matrix *mA, *mB, *mC;
219
type = pdp_packet_matrix_get_type(pC);
220
mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type);
221
mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB,type);
222
mC = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pC,type);
224
if (!(mA && mB)) return 1;
228
case PDP_MATRIX_TYPE_RFLOAT:
229
return gsl_blas_sgemm(TransA, TransB, scale_r, (gsl_matrix_float *)mA,
230
(gsl_matrix_float *)mB, 1.0f, (gsl_matrix_float *)mC);
231
case PDP_MATRIX_TYPE_RDOUBLE:
232
return gsl_blas_dgemm(TransA, TransB, (double)scale_r, (gsl_matrix *)mA,
233
(gsl_matrix *)mB, 1.0, (gsl_matrix *)mC);
234
case PDP_MATRIX_TYPE_CFLOAT:
235
return gsl_blas_cgemm(TransA, TransB, cf_scale, (gsl_matrix_complex_float *)mA,
236
(gsl_matrix_complex_float *)mB, cf_one, (gsl_matrix_complex_float *)mC);
237
case PDP_MATRIX_TYPE_CDOUBLE:
238
return gsl_blas_zgemm(TransA, TransB, cd_scale, (gsl_matrix_complex *)mA,
239
(gsl_matrix_complex *)mB, cd_one, (gsl_matrix_complex *)mC);
245
/* matrix vector mul: C is defining type
246
returns 0 on success */
247
int pdp_packet_matrix_blas_mv
249
CBLAS_TRANSPOSE_t TransA,
257
gsl_complex_float cf_scale = {{scale_r, scale_i}};
258
gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}};
259
gsl_complex_float cf_one = {{1.0f, 0.0f}};
260
gsl_complex cd_one = {{1.0, 0.0}};
264
type = pdp_packet_matrix_get_type(pA);
265
mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type);
266
vb = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pb,type);
267
vc = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pc,type);
269
if (!(vb && vc)) return 1;
273
case PDP_MATRIX_TYPE_RFLOAT:
274
return gsl_blas_sgemv(TransA, scale_r, (gsl_matrix_float *)mA,
275
(gsl_vector_float *)vb, 1.0f, (gsl_vector_float *)vc);
276
case PDP_MATRIX_TYPE_RDOUBLE:
277
return gsl_blas_dgemv(TransA, (double)scale_r, (gsl_matrix *)mA,
278
(gsl_vector *)vb, 1.0, (gsl_vector *)vc);
279
case PDP_MATRIX_TYPE_CFLOAT:
280
return gsl_blas_cgemv(TransA, cf_scale, (gsl_matrix_complex_float *)mA,
281
(gsl_vector_complex_float *)vb, cf_one, (gsl_vector_complex_float *)vc);
282
case PDP_MATRIX_TYPE_CDOUBLE:
283
return gsl_blas_zgemv(TransA, cd_scale, (gsl_matrix_complex *)mA,
284
(gsl_vector_complex *)vb, cd_one, (gsl_vector_complex *)vc);
292
t_pdp_symbol *_pdp_matrix_get_description(t_pdp *header)
294
t_matrix *m = (t_matrix *)(&header->info.raw);
295
char description[100];
296
char *c = description;
299
if (!header) return pdp_gensym("invalid");
300
else if (!header->desc){
301
/* if description is not defined, try to reconstruct it (for backwards compat) */
302
if (header->type == PDP_MATRIX){
304
c += sprintf(c, "matrix");
307
case PDP_MATRIX_TYPE_RFLOAT: c += sprintf(c, "/float/real"); break;
308
case PDP_MATRIX_TYPE_CFLOAT: c += sprintf(c, "/float/complex"); break;
309
case PDP_MATRIX_TYPE_RDOUBLE: c += sprintf(c, "/double/real"); break;
310
case PDP_MATRIX_TYPE_CDOUBLE: c += sprintf(c, "/double/complex"); break;
312
c += sprintf(c, "/unknown"); goto exit;
315
c += sprintf(c, "/%dx%d", (int)m->rows, (int)m->columns);
318
return pdp_gensym(description);
320
else return pdp_gensym("unknown");
322
else return header->desc;
327
static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src);
328
static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src);
329
static void _pdp_matrix_reinit(t_pdp *dst);
330
static void _pdp_matrix_cleanup(t_pdp *dst);
332
static size_t _pdp_matrix_mdata_byte_size(u32 rows, u32 columns, u32 type)
334
size_t dsize = rows * columns;
336
case PDP_MATRIX_TYPE_RFLOAT: dsize *= sizeof(float); break;
337
case PDP_MATRIX_TYPE_CFLOAT: dsize *= 2*sizeof(float); break;
338
case PDP_MATRIX_TYPE_RDOUBLE: dsize *= sizeof(double); break;
339
case PDP_MATRIX_TYPE_CDOUBLE: dsize *= 2*sizeof(double); break;
346
static size_t _pdp_matrix_pdata_vector_size(u32 rows, u32 columns)
348
return (rows>columns ? rows : columns);
352
static void _pdp_matrix_init(t_pdp *dst, u32 rows, u32 columns, u32 type)
355
t_matrix *m = (t_matrix *)(&dst->info.raw);
356
void *d = ((void *)dst) + PDP_HEADER_SIZE;
357
int matrix_bytesize = _pdp_matrix_mdata_byte_size(rows, columns, type);
358
int permsize = _pdp_matrix_pdata_vector_size(rows, columns);
363
m->columns = columns;
365
/* set the block data */
366
m->block.size = matrix_bytesize;
367
m->block.data = (double *)d;
369
/* set the vector view */
370
m->vector.size = (1==columns) ? rows : columns;
371
m->vector.stride = 1;
372
m->vector.data = (double *)d;
373
m->vector.block = &m->block;
376
/* set the matrix view */
377
m->matrix.gsl_rows = rows;
378
m->matrix.gsl_columns = columns;
379
m->matrix.tda = columns;
380
m->matrix.data = (double *)d;
381
m->matrix.block = &m->block;
384
/* set the permutation object & init */
385
m->perm.size = permsize;
386
m->perm.data = (size_t *)(d + matrix_bytesize);
387
for(i=0; i<permsize; i++) m->perm.data[i] = i;
390
/* init packet header */
391
dst->theclass = matrix_class;
392
dst->desc = _pdp_matrix_get_description(dst);
396
static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src)
398
t_matrix *m = (t_matrix *)(&src->info.raw);
399
_pdp_matrix_init(dst, m->rows, m->columns, m->type);
402
static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src)
404
_pdp_matrix_clone(dst, src);
405
memcpy(dst + PDP_HEADER_SIZE, src + PDP_HEADER_SIZE, src->size - PDP_HEADER_SIZE);
406
//post("matrix copy successful");
408
static void _pdp_matrix_reinit(t_pdp *dst)
410
/* nothing to do, assuming data is correct */
412
static void _pdp_matrix_cleanup(t_pdp *dst)
414
/* no extra memory to free */
417
int pdp_packet_new_matrix(u32 rows, u32 columns, u32 type)
422
/* compute the blocksize for the matrix data */
423
/* if 0, something went wrong -> return invalid packet */
424
if (!(dsize = _pdp_matrix_mdata_byte_size(rows, columns, type))) return -1;
425
dsize += sizeof(size_t) * _pdp_matrix_pdata_vector_size(rows, columns);
427
p = pdp_packet_new(PDP_MATRIX, dsize);
428
if (-1 == p) return -1;
429
header = pdp_packet_header(p);
431
_pdp_matrix_init(header, rows, columns, type);
436
int pdp_packet_new_matrix_product_result(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, int pA, int pB)
439
u32 type = pdp_packet_matrix_get_type(pA);
441
u32 colA, colB, rowA, rowB;
443
/* check if A is a matrix */
444
if (!type) return -1;
446
mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA, type);
447
mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB, type);
449
/* check if A and B are same type */
453
if (TransA == CblasNoTrans){
455
colA = mA->gsl_columns;
458
rowA = mA->gsl_columns;
463
if (TransB == CblasNoTrans){
465
colB = mB->gsl_columns;
468
rowB = mB->gsl_columns;
472
/* check if sizes are compatible */
473
if (colA != rowB) return -1;
475
/* create new packet */
476
return pdp_packet_new_matrix(rowA, colB, type);
479
void pdp_packet_matrix_setzero(int p)
482
if (!pdp_packet_matrix_isvalid(p)) return;
483
m = (t_matrix *) pdp_packet_subheader(p);
484
memset(m->block.data, 0, m->block.size);
489
/* type conversion programs */
490
static int _pdp_packet_convert_matrix_to_greyimage(int packet, t_pdp_symbol *dest_template)
492
t_pdp *header = pdp_packet_header(packet);
493
t_matrix *matrix = (t_matrix *)pdp_packet_subheader(packet);
494
void *data = pdp_packet_data(packet);
504
new_p = pdp_packet_new_image_grey(c,r);
505
if (-1 == new_p) return -1;
506
new_data = (s16 *)pdp_packet_data(new_p);
508
/* convert first channel */
509
switch(matrix->type){
510
case PDP_MATRIX_TYPE_RFLOAT:
511
for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i] * (float)0x8000); break;
512
case PDP_MATRIX_TYPE_CFLOAT: //only copy real channel
513
for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i<<1] * (float)0x8000); break;
514
case PDP_MATRIX_TYPE_RDOUBLE:
515
for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i] * (float)0x8000); break;
516
case PDP_MATRIX_TYPE_CDOUBLE: //only copy real channel
517
for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i<<1] * (float)0x8000); break;
519
pdp_packet_mark_unused(new_p);
525
static int _pdp_packet_convert_image_to_matrix(int packet, t_pdp_symbol *dest_template, int type)
527
t_pdp *header = pdp_packet_header(packet);
528
t_image *image = pdp_packet_image_info(packet);
529
s16 *data = (s16 *)pdp_packet_data(packet);
533
int encoding = image->encoding;
536
if (!pdp_packet_image_isvalid(packet)) return -1;
541
new_p = pdp_packet_new_matrix(h,w, type);
542
if (-1 == new_p) return -1;
543
new_data = pdp_packet_data(new_p);
549
/* convert first channel */
551
case PDP_MATRIX_TYPE_RFLOAT:
552
for (i=0; i<nbelements; i++)
553
((float *)new_data)[i] = ((float)data[i]) * (1.0f / (float)0x8000); break;
554
case PDP_MATRIX_TYPE_RDOUBLE:
555
for (i=0; i<nbelements; i++)
556
((double *)new_data)[i] = ((double)data[i]) * (1.0f / (double)0x8000); break;
557
case PDP_MATRIX_TYPE_CFLOAT:
558
for (i=0; i<nbelements; i++){
559
((float *)new_data)[i*2] = ((float)data[i]) * (1.0f / (float)0x8000);
560
((float *)new_data)[i*2+1] = 0.0f;
563
case PDP_MATRIX_TYPE_CDOUBLE:
564
for (i=0; i<nbelements; i++){
565
((double *)new_data)[i*2] = ((double)data[i]) * (1.0f / (double)0x8000);
566
((double *)new_data)[i*2+1] = 0.0;
570
pdp_post("_pdp_packet_convert_image_to_matrix: INTERNAL ERROR");
574
pdp_packet_mark_unused(new_p);
583
static int _pdp_packet_convert_image_to_floatmatrix(int packet, t_pdp_symbol *dest_template){
584
return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RFLOAT);}
585
static int _pdp_packet_convert_image_to_doublematrix(int packet, t_pdp_symbol *dest_template){
586
return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RDOUBLE);}
587
static int _pdp_packet_convert_image_to_complexfloatmatrix(int packet, t_pdp_symbol *dest_template){
588
return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CFLOAT);}
589
static int _pdp_packet_convert_image_to_complexdoublematrix(int packet, t_pdp_symbol *dest_template){
590
return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CDOUBLE);}
593
static int _pdp_packet_matrix_convert_fallback(int packet, t_pdp_symbol *dest_template)
595
pdp_post("can't convert image type %s to %s",
596
pdp_packet_get_description(packet)->s_name, dest_template->s_name);
602
/* this seems like a nice spot to place a dummy gsl signal handler */
603
static void _gsl_error_handler (const char * reason,
608
//pdp_post("gsl error:\nREASON: %s\nFILE:%s\nLINE:%d\nERRNO:%d", reason, file, line, gsl_errno);
611
static int pdp_matrix_factory(t_pdp_symbol *type)
615
void pdp_matrix_setup(void)
617
t_pdp_conversion_program *program;
620
/* setup the class */
621
matrix_class = pdp_class_new(pdp_gensym("matrix/*/*/*"), pdp_matrix_factory);
622
matrix_class->copy = _pdp_matrix_copy;
623
//matrix_class->clone = _pdp_matrix_clone; // is now solved through (symbol based) constructor
624
matrix_class->wakeup = _pdp_matrix_reinit;
625
matrix_class->cleanup = _pdp_matrix_cleanup;
628
/* image -> matrix: default = double/real (most ops available) */
629
program = pdp_conversion_program_new(_pdp_packet_convert_image_to_doublematrix, 0);
630
pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/real/*"), program);
631
program = pdp_conversion_program_new(_pdp_packet_convert_image_to_floatmatrix, 0);
632
pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/real/*"), program);
633
program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexdoublematrix, 0);
634
pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/complex/*"), program);
635
program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexfloatmatrix, 0);
636
pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/complex/*"), program);
638
/* matrix -> image */
639
program = pdp_conversion_program_new(_pdp_packet_convert_matrix_to_greyimage, 0);
640
pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/grey/*"), program);
643
program = pdp_conversion_program_new(_pdp_packet_matrix_convert_fallback, 0);
644
pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/*/*"), program);
645
pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("bitmap/*/*"), program);
646
pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/*/*/*/*"), program);
647
pdp_type_register_conversion(pdp_gensym("bitmap/*/*"), pdp_gensym("matrix/*/*/*/*"), program);
649
/* setup gsl handler */
650
if(gsl_set_error_handler(_gsl_error_handler)){
651
pdp_post("pdp_matrix_setup: WARNING: overriding gsl error handler.");
added
removed
system/type/pdp_matrix.c
