Inversa de una matriz e intersección de 2 planos

[Findeton]

                                Hola,estoy intentando crear una serie de funciones que me sirvan para detectar la colisión entre polígonos, esferas y líneas y me he encontrado un problema:

Estoy intentando encontrar la intersección entre 2 planos, los cuales han sido calculados a partir de 2 polígonos, dándome así la ecuación de cada plano Ax +By + Cz +D = 0, siendo el vector normal de cada plano (x=A, y=B, z=C). Entonces, utilizando el método de Microsoft, que viene descrito en la web http://research.microsoft.com/~jckrumm/Int..._two_planes.htm , sólo necesito resolver un sistema de ecuaciones en forma de matrices.


Para ello, me fui a la web http://www.sc.ehu.es/sbweb/fisica/cursoJav...iz1/matriz1.htm donde explican que para saber la solución de A*S=B (donde A es una matriz conocida, B un vector conocido y S el vector solución), hay que calcular la inversa de A y multiplicarla por B, y te dará S. Entonces, como en esa misma web me viene el código para calcular la inversa... pensé que todo estaba solucionado... pasé el código a C++ y... pues sólo calcula bien la inversa de matrices 2x2, cuando el método descrito en la web es para nxn, donde n es cualquier número entero.

Así que posteo el código para que comprobéis si hay algún error, si podéis  :(.

Ésta es la función que calcula la matriz inversa:

Código [Seleccionar] Expand



Matriz matriz_inversa(Matriz d){

       int n=d.n;     //dimensión de la matriz

 Matriz a(n);

 a=d;

       Matriz b(n);    //matriz de los términos independientes

       Matriz c(n);    //matriz de las incógnitas



//matriz unidad

       for(int i=0; i<n; i++){

           b.x[i][i]=1.0;

       }

//transformación de la matriz y de los términos independientes

       for(int k=0; k<n-1; k++){

           for(int i=k+1; i<n; i++){

//términos independientes

               for(int s=0; s<n; s++){

                   b.x[i][s]-=a.x[i][k]*b.x[k][s]/a.x[k][k];

               }

//elementos de la matriz

               for(int j=k+1; j<n; j++){

                   a.x[i][j]-=a.x[i][k]*a.x[k][j]/a.x[k][k];

               }

           }

       }

//cálculo de las incógnitas, elementos de la matriz inversa

       for(int s=0; s<n; s++){

           c.x[n-1][s]=b.x[n-1][s]/a.x[n-1][n-1];

           for(int i=n-2; i>=0; i--){

               c.x[i][s]=b.x[i][s]/a.x[i][i];

               for(int k=n-1; k>i; k--){

                   c.x[i][s]-=a.x[i][k]*c.x[k][s]/a.x[i][i];

               }

           }

       }

       return c;

   }



Y aquí está la definicíon de la clase Matriz:

Código [Seleccionar] Expand



class Matriz {



public:

int n;         //dimensión

   double** x;     //el valor

   Matriz(int m); //constructor

  ~Matriz();     //destructor

  void mostrar();

   };





//Aquí definimos el constructor



Matriz::Matriz(int m) {

       n=m;

       x= new double*[n];

       for(int i=0; i<n; i++){

  x[i]= new double[n];

 }  

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

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

   }



//Aquí definimos el destructor



Matriz::~Matriz() {

   

                     

   }



void Matriz::identidad(){

       for(int i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

 for(i=0; i<n; i++)

  x[i][i]=1;

}

void Matriz::identidad(double m){

       for(int i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

 for(i=0; i<n; i++)

  x[i][i]=m;

}



void Matriz::mostrar(){



for(int i=0; i<n;i++){ //hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

 for(int k=0; k<n;k++){//hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

        cout << x[i][k] << " ";

 }

 cout << endl;

}

}






Si necesitáis más datos, por favor pedirlos, pero ya os digo que un código de ejemplo sería:

Código [Seleccionar] Expand



#include "iostream.h"



/*Poned aquí la definición de la matriz y sus funciones*/



int main(int argc, char* argv[])

{

Matriz b(3);



b.x[0][0]=4;

b.x[0][1]=1;

b.x[0][2]=1;

b.x[1][0]=2;

b.x[1][1]=3;

b.x[1][2]=1;

b.x[2][0]=2;

b.x[2][1]=3;

b.x[2][2]=1;



   

b=matriz_inversa(b);

b.mostrar();



return 0;

}