1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
/*
* iemmatrix
*
* objects for manipulating simple matrices
* mostly refering to matlab/octave matrix functions
*
* Copyright (c) IOhannes m zmölnig, forum::für::umläute
* IEM, Graz, Austria
*
* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
*
*/
#include "iemmatrix.h"
/* mtx_inverse */
static t_class *mtx_inverse_class;
t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol, int*error){
/*
* row==col==rowclo
* input=t_matrixfloat[row*col]
* output=t_matrixfloat[row*col]
*/
int i, k;
t_matrixfloat *a1, *b1, *a2, *b2;
int ok=0; /* error counter */
int col=rowcol, row=rowcol, row2=row*col;
t_matrixfloat *original=input;
t_matrixfloat *inverted = 0;
if(input==0)return 0;
/* 1a reserve space for the inverted matrix */
inverted=(t_matrixfloat *)getbytes(sizeof(t_matrixfloat)*row2);
if(inverted==0)return 0;
/* 1b make an eye-shaped float-buf for B */
i=row2;
b1=inverted;
while(i--)*b1++=0;
i=row;
b1=inverted;
while(i--)b1[i*(row+1)]=1;
/* 2. do the Gauss-Jordan */
for (k=0;k<row;k++) {
/* adjust current row */
t_matrixfloat diagel = original[k*(col+1)];
t_matrixfloat i_diagel = diagel?1./diagel:0;
if (!diagel)ok++;
/* normalize current row (set the diagonal-element to 1 */
a2=original+k*col;
b2=inverted+k*col;
i=row;
while(i--){
*a2++*=i_diagel;
*b2++*=i_diagel;
}
/* eliminate the k-th element in each row by adding the weighted normalized row */
a2=original+k*row;
b2=inverted+k*row;
for(i=0;i<row;i++)
if (i-k) {
t_matrixfloat f=-*(original+i*row+k);
int j = row;
a1=original+i*row;
b1=inverted+i*row;
while (j--) {
*(a1+j)+=f**(a2+j);
*(b1+j)+=f**(b2+j);
}
}
}
if (ok)post("mtx_inverse: couldn't really invert the matrix !!! %d error%c", ok, (ok-1)?'s':0);
if(error!=0)*error=ok;
return inverted;
}
static void mtx_inverse_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
{
/* maybe we should do this in double or long double ? */
int row=atom_getfloat(argv);
int col=atom_getfloat(argv+1);
int error=0;
t_matrixfloat *original, *inverted;
if(row*col+2>argc){
post("mtx_print : sparse matrices not yet supported : use \"mtx_check\"");
return;
}
/* reserve memory for outputting afterwards */
adjustsize(x, col, row);
/* 1. extract values of A to float-buf */
original=matrix2float(argv);
if (row==col){
/* fine, the matrix is square */
inverted=mtx_doInvert(original, row, &error);
} else {
/* we'll have to do the pseudo-inverse:
* P=A'*inv(A*A') if row<col
* P=inv(A'*A)*A' if col<row
*/
t_matrixfloat*transposed, *invertee;
int inverteeCol=0;
transposed=mtx_doTranspose(original, row, col);
if(row>col){
inverteeCol=col;
invertee =mtx_doMultiply(col, transposed, row, original, col);
inverted =mtx_doMultiply(col, mtx_doInvert(invertee, col, &error), col, transposed, row);
} else {
inverteeCol=row;
invertee =mtx_doMultiply(row, original, col, transposed, row);
inverted =mtx_doMultiply(col, transposed, row, mtx_doInvert(invertee, row, &error), row);
}
freebytes(transposed, sizeof(t_matrixfloat)*col*row);
freebytes(invertee , sizeof(t_matrixfloat)*inverteeCol*inverteeCol);
}
/* 3. output the matrix */
/* 3a convert the floatbuf to an atombuf; */
float2matrix(x->atombuffer, inverted);
/* 3b destroy the buffers */
freebytes(original, sizeof(t_matrixfloat)*row*col);
if(error)outlet_bang(x->x_outlet);
/* 3c output the atombuf; */
matrix_bang(x);
}
static void *mtx_inverse_new(t_symbol *s, int argc, t_atom *argv)
{
t_matrix *x = (t_matrix *)pd_new(mtx_inverse_class);
outlet_new(&x->x_obj, 0);
x->col=x->row=0;
x->atombuffer=0;
x->x_outlet=outlet_new(&x->x_obj, 0);
return (x);
}
void mtx_inverse_setup(void)
{
mtx_inverse_class = class_new(gensym("mtx_inverse"), (t_newmethod)mtx_inverse_new,
(t_method)matrix_free, sizeof(t_matrix), 0, A_GIMME, 0);
class_addbang (mtx_inverse_class, matrix_bang);
class_addmethod(mtx_inverse_class, (t_method)mtx_inverse_matrix, gensym("matrix"), A_GIMME, 0);
class_sethelpsymbol(mtx_inverse_class, gensym("iemmatrix/mtx_inverse"));
}
void iemtx_inverse_setup(void){
mtx_inverse_setup();
}
|
