[C] matrice inversa 4x4 con algoritmo di Gauss

[Thunderx]

ragazzi sto veramente uscendo pazzo......
voglio implementare l'algoritmo di gauss per calcolare la matrice inversa di una 4x4, ma non so più che pesci prendere visto che i risultati non sono mai corretti...il listato al post 3(Ogni consiglio è ben accetto)

[Thunderx]

edit

[Thunderx]

#include <math.h>
#include <stdio.h>
void inversa(float (*a)[4]) {

	int i;
	int j;
        int z;
	float b[4][4] = {{ 1, 0, 0, 0},{ 0, 1, 0, 0},{ 0, 0, 1, 0},{ 0, 0, 0, 1}};
	float c;
        float d;
        for (z=0;z<3;z++) {
	for (i=z;i<3;i++) {
		if (a[i+1][z]=0){
                     continue;}
else
             {    c = -(a[i+1][z] / a[i][z]);
                for (j=z;j<3;j++) {
                a[i + 1][j] = (c * (a[i][j])) + a[i + 1][j];
		b[i + 1][j] = (c * (b[i][j])) + b[i + 1][j];
}
}	
}
}
  for (z=0;z<3;z++) {
	for (i=(3-z);i>0;i--) {
if(a[i-1][3-z]=0){
continue; }

else
{
d = -(a[i-1][(3-z)]) / a[i][3-z];
	               for (j=(3-z);j>0;j--)
 {		
                        a[i-1][j]=(d*(a[i][j]))+a[i-1][j];
			b[i-1][j]=(d*(b[i][j]))+b[i-1][j];
}
}	
}
}
	for (i = 0; i < 4; i++) {
		for (j = 0; j < 4; j++) {
			b[i][j] = b[i][j] / a[i][i];
			a[i][j] = a[i][j] / a[i][i];
		}
	}
	for (i = 0; i <4; i++) {
		printf("\n");
		for (j = 0; j <4; j++) {
			printf("%f" , a[i][j]);
		}
	}
	for (i = 0; i < 4; i++) {
		printf("\n");
		for (j = 0; j < 4; j++) {
                      printf("%f", b[i][j]);
			}
		}
}


int main(void) {
		float mat[4][4];
		printf("Scrivi gli elementi della matrice:\n");
		int i;
		int j;
            for (i = 0;i < 4; i++) {
			for (j = 0; j < 4 ;j ++) {
				scanf("%f", &mat[i][j]);
                     }
		}
inversa(mat);
return (0);
}

ecco l'uscita:

1.000000nannannan
nannannannan
nannannannan
nannannan1.000000
1.000000nannannan
nannannannan
nannannannan

[altairz]

...
if (a[i+1][z] = 0)
...
if (a[i-1][3-z] = 0)
...

credo che un bug sia in questo punto

[Thunderx]

Originariamente inviato da altairz

...
if (a[i+1][z] = 0)
...
if (a[i-1][3-z] = 0)
...

credo che un bug sia in questo punto

alla fine quel pezzo l'ho tolto modificando la condizione sul coefficiente.....ecco il listato della parte modificata

       for (z=0;z<3;z++) {
	for (i=z;i<3;i++) {
   c = -(a[i+1][z] / a[z][z]);
                for (j=z;j<4;j++) {
                 
                a[i + 1][j] = c*a[i][j] + a[i + 1][j];
		
                b[i + 1][j] =c*b[i][j]+ b[i + 1][j];
}
}	
}
  for (z=0;z<3;z++) {
	for (i=(3-z);i>0;i--){
                   d = -(a[i-1][(3-z)]) / a[3-z][3-z];
	               for (j=(3-z);j> -1;j--){
                        a[i-1][j]=d*a[i][j]+a[i-1][j];
                        b[i-1][j]=d*b[i][j]+b[i-1][j];
}	
}
}

questa è l'uscita

1.000000-1.0000002.000000-1.000000
0.0000001.000000-1.0000001.000000
0.0000000.0000001.000000-1.000000
0.0000000.0000000.0000001.000000

con ingresso

1 1 1 1
0 1 1 1
0 0 1 1
0 0 0 1

[Thunderx]

mi sono accorto di un piccolo errore nell'implementazione dell'algoritmo e ho nuovamente modificato il programma

      for (z=0;z<3;z++) {
	for (i=z;i<3;i++) {
   c = -(a[i+1][z] / a[z][z]);
                for (j=z;j<4;j++) {
                a[i][j]= c * a[i][j]; 
                a[i + 1][j] = a[i][j] + a[i + 1][j];
		b[i][j]= c * b[i][j];
                b[i + 1][j] =b[i][j]+ b[i + 1][j];
}}}
  for (z=0;z<3;z++) {
	for (i=(3-z);i>0;i--){
                   d = -(a[i-1][(3-z)]) / a[3-z][3-z];
	               for (j=(3-z);j> -1;j--)
 {		        a[i][j]= d*a[i][j];
                        a[i-1][j]=a[i][j]+a[i-1][j];
			b[i][j]= d*b[i][j];
                        b[i-1][j]=b[i][j]+b[i-1][j];

}}}

con lo stesso ingresso l'uscita è

nannannannan
nannannannan
nannannannan
nannannannan
nannannannan
nannannannan
nannannannan
nannannannan

come mai?eppure dovrebbe essere giusto

[UnknownSoldier]

Perchè non usi il tag [code] o [php] ? Renderebbe la lettura del codice molto più leggibile.

[Thunderx]

Originariamente inviato da UnknownSoldier

Perchè non usi il tag [code] o [php] ? Renderebbe la lettura del codice molto più leggibile.

ok metto subito a posto. grazie per la segnalazione

[sottovento]

Se ottieni NaN e' perche' vai a fare le divisioni per zero.
In effetti hai un paio di divisioni che andrebbero controllate. Cosa devi fare, secondo l'algoritmo, nel caso che il divisore sia zero?

[Thunderx]

Originariamente inviato da sottovento

Se ottieni NaN e' perche' vai a fare le divisioni per zero.
In effetti hai un paio di divisioni che andrebbero controllate. Cosa devi fare, secondo l'algoritmo, nel caso che il divisore sia zero?

in teoria non dovrebbe mai verificarsi una divisione per 0... anche perchè il divisore è sempre il pivot che per definizione è non nullo

[sottovento]

Questo vuol dire che puoi aggiungere una assert al tuo codice per essere veramente sicuro di questo

[||ElChE||88]

Originariamente inviato da sottovento

Se ottieni NaN e' perche' vai a fare le divisioni per zero.

Se non sbaglio dividendo per 0 un numero a virgola mobile (IEEE 754) il risultato è +/-∞, non NaN.

[Thunderx]

Originariamente inviato da sottovento

Questo vuol dire che puoi aggiungere una assert al tuo codice per essere veramente sicuro di questo

che cos'è una assert?

[sottovento]

si scrive assert(condizione che vuoi verificare) e se la condizione non e' verificata, si arrabbia e fa terminare il programma.

Hai provato a far girare il debugger?

[sottovento]

Originariamente inviato da ||ElChE||88

Se non sbaglio dividendo per 0 un numero a virgola mobile (IEEE 754) il risultato è +/-∞, non NaN.

Penso che tu abbia ragione, ma non sono sicuro al 100%.
Spero che questo esercizio me lo chiarisca

[Thunderx]

Originariamente inviato da sottovento

si scrive assert(condizione che vuoi verificare) e se la condizione non e' verificata, si arrabbia e fa terminare il programma.

Hai provato a far girare il debugger?

ho fatto una nuov revision del codice

		}     
        for (z=0;z<3;z++) {
	for (i=z;i<3;i++) {
  if (z=0){ c = -(a[i+1][z] / a[z][z]);}   
else {c=-(a[i+1][z] / p[z][z]);}
                for (j=z;j<4;j++) {
               p[i][j]=c*a[i][j]; 
               a[i + 1][j] =p[i][j] + a[i + 1][j];
	       q[i][j]=c*b[i][j];
               b[i + 1][j] =q[i][j] + b[i + 1][j];
}
}	
}
  for (z=0;z<3;z++) {
	for (i=(3-z);i>0;i--){
                if(z=0) {   d = -(a[i-1][(3-z)]) / a[3-z][3-z];}
else{d = -(a[i-1][(3-z)]) / p[3-z][3-z];}
	               for (j=(3-z);j> -1;j--){
                       p[i][j]=d*a[i][j];
                       a[i-1][j]=p[i][j]+a[i-1][j];
                       q[i][j]=d*b[i][j];
                       b[i-1][j]=q[i][j]+b[i-1][j];
}	
}
}
	for (i = 0; i < 4; i++) {
		for (j = 0; j < 4; j++) {
			b[i][j] = b[i][j] / a[i][i];
			a[i][j] = a[i][j] / a[i][i];
		}
	}

il mio problema era quello di aggiornare i valori delle matrice a , moltiplicandola per i vari coefficienti di volta in volta e quindi ho inserito l'if

[sottovento]

Quell'if non va bene (z=0 e' un assegnamento).

Ho fatto un giro di debugger fermandomi immediatamente al primo ciclo (il risultato e' gia' incorretto li', per cui non ho proseguito).

Dopo la prima iterazione del ciclo

        for (i=z;i<3;i++)
        {

il valore di a[0][0] e' zero! Questo e' il pivot per il ciclo successivo (in cui z e' ancora uguale a zero). Quindi vai poi a dividere per questo valore.

Wikipedia suggerisce una cosa cosi':

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k
    end if
  end for
  if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]
    Now A[i,j] will have the value 1.
    for u := i+1 to m do
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1
  end if
  j := j + 1
end while

[Thunderx]

Originariamente inviato da sottovento

Quell'if non va bene (z=0 e' un assegnamento).

Ho fatto un giro di debugger fermandomi immediatamente al primo ciclo (il risultato e' gia' incorretto li', per cui non ho proseguito).

Dopo la prima iterazione del ciclo

        for (i=z;i<3;i++)
        {

il valore di a[0][0] e' zero! Questo e' il pivot per il ciclo successivo (in cui z e' ancora uguale a zero). Quindi vai poi a dividere per questo valore.

Wikipedia suggerisce una cosa cosi':

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k
    end if
  end for
  if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]
    Now A[i,j] will have the value 1.
    for u := i+1 to m do
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1
  end if
  j := j + 1
end while

scusami forse non mi sono spiegato,a[i][j] è una matrice che viene definita dall'utente nel main.Visto che poi dovrà essere collegata ad un programma più grande dove la matrice a difficilmente avrà termini nulli.
la matrice di prova che io uso è
1111
0111
0011
0001

[sottovento]

Originariamente inviato da Thunderx

scusami forse non mi sono spiegato,a[i][j] è una matrice che viene definita dall'utente nel main.Visto che poi dovrà essere collegata ad un programma più grande dove la matrice a difficilmente avrà termini nulli.
la matrice di prova che io uso è
1111
0111
0011
0001

Si, l'avevo capito ed ho usato la stessa matrice per fare un po' di debug. Ma poi, ciclo dopo ciclo, vai a modificare quella matrice e vai a scegliere un valore che e' uguale a zero per effettuare la divisione.

Ho completamente scordato questo algoritmo, per cui perdona la mia domanda: sei sicuro che devi modificare la matrice "on-the-fly" (i.e. durante l'esecuzione del ciclo) e non alla fine?
Normalmente questo tipo di algoritmi va a preparare una seconda matrice on i risultati. Questa matrice poi diventa quella di partenza per il ciclo successivo...

Inoltre, per quel poco che ho letto su wikipedia, sembra che sia normale il fatto che il tuo elemento sia azzerato, proprio per come e' fatto l'algoritmo. Forse e' semplicemente la scelta del pivot che va cambiata.

Spero ti sia di aiuto....

[Thunderx]

Originariamente inviato da sottovento

Si, l'avevo capito ed ho usato la stessa matrice per fare un po' di debug. Ma poi, ciclo dopo ciclo, vai a modificare quella matrice e vai a scegliere un valore che e' uguale a zero per effettuare la divisione.

Ho completamente scordato questo algoritmo, per cui perdona la mia domanda: sei sicuro che devi modificare la matrice "on-the-fly" (i.e. durante l'esecuzione del ciclo) e non alla fine?
Normalmente questo tipo di algoritmi va a preparare una seconda matrice on i risultati. Questa matrice poi diventa quella di partenza per il ciclo successivo...

Inoltre, per quel poco che ho letto su wikipedia, sembra che sia normale il fatto che il tuo elemento sia azzerato, proprio per come e' fatto l'algoritmo. Forse e' semplicemente la scelta del pivot che va cambiata.

Spero ti sia di aiuto....

preso da wikipedia

L'algoritmo di Gauss trasforma una qualsiasi matrice in una matrice a scalini tramite mosse di Gauss. Funziona nel modo seguente:

1. Se la prima riga ha il primo elemento nullo, scambiala con una riga che ha il primo elemento non nullo. Se tutte le righe hanno il primo elemento nullo, vai al punto 3.
2. Per ogni riga Ai con primo elemento non nullo, eccetto la prima (i > 1), moltiplica la prima riga per un coefficiente scelto in maniera tale che la somma tra la prima riga e Ai abbia il primo elemento nullo (quindi coefficiente − Ai1 / A11). Sostituisci Ai con la somma appena ricavata.
3. Adesso sulla prima colonna tutte le cifre, eccetto forse la prima, sono nulle. A questo punto ritorna al punto 1 considerando la sottomatrice che ottieni cancellando la prima riga e la prima colonna. Le successive mosse di Gauss (scambi, moltiplicazioni e somme) andranno comunque fatte su tutta la matrice.

quindi come vedi l'algoritmo modifica il coefficiente c per ogni incremento di i e naturalmente anche di z.
adesso con l'ultimo codice scritto non ho più problemi di nan, ma solodi 3 valori su 16 sballatidovuti al fatto che la riga z-esima (la prima di ogni sottomatrice dell'algoritmo) non viene mai moltiplicata per c ma rimane statica
Scusa ma perchè il ciclo if di prima non andrebbe bene?come dovrebbe essere?
ti volevo ringraziare in ogni caso, perchè a occhio dalle tue parti dovrebbe essere parecchio tardi