max/msp gen~ 数値計算　逆行列を求める。

max/msp gen~でガウスジョルダン法を使って逆行列を求めます。

このような感じで3行3列の行列の要素を入力すると逆行列を出力します。

MATLABで同じ行列の逆行列を求めたものです。
おおよそ同じような解が出力されています。

gen~の中身はこんな感じです。

gen~の中身です。

Data a(3,3);
Data b(3,6);

N=3;eps=0.00001;
piv0t=0;delta=0;x=0;y=0;i=0;j=0;k=0;
no_exist=0;
b_ij=0;b_kj;

poke(a,in1,0,0);poke(a,in2,0,1);poke(a,in3,0,2);
poke(a,in4,1,0);poke(a,in5,1,1);poke(a,in6,1,2);
poke(a,in7,2,0);poke(a,in8,2,1);poke(a,in9,2,2);

for(y=0;y<N;y+=1){
  for(x=0;x<N;x+=1){
    poke(b,peek(a,y,x),y,x);
    poke(b,0.0,y,x+N);
  }
  poke(b,1.0,y,y+N);	
}

for(i=0;i<N;i+=1){
 pivot = peek(b,i,i);
 if( abs(pivot) < eps){
      no_exist=1;
   }	
 
 for(j=1; j<2*N;j+=1){
   b_ij=peek(b,i,j);
   poke(b,b_ij/pivot,i,j);
 }
 for(k=0;k<N;k+=1){
   delta = peek(b,k,i);
   for(j=i;j<2*N;j+=1){
     if(k != i){
      b_kj=peek(b,k,j)-delta*peek(b,i,j);
      poke(b,b_kj,k,j);
     }
   }
 }
}

out1=peek(b,0,0+N);out2=peek(b,0,1+N);out3=peek(b,0,2+N);
out4=peek(b,1,0+N);out5=peek(b,1,1+N);out6=peek(b,1,2+N);
out7=peek(b,2,0+N);out8=peek(b,2,1+N);out9=peek(b,2,2+N); out10=no_exist;

4 ガウス・ジョルダン法

ガウスの消去法 - Wikipedia