[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):D

[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]

...
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

:muro: :muro: :muro: :muro:
come mai?eppure dovrebbe essere giusto

[UnknownSoldier]

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

[Thunderx]

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]

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 :D

[||ElChE||88]

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]

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

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]

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]

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]

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]

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]

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

[sottovento]

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 e naturalmente anche di z.
Scusa ma perchè il ciclo if di prima non andrebbe bene?

Hai ragione, ma probabilmente non ci siamo capiti.
Quando devi calcolare c, utilizzi un elemento sulla diagonale per fare la divisione, confidente del fatto che la tua matrice di test ha la diagonale non nulla:

c = -(a[i+1][z] / a[z][z]);

In realta', questo e' vero per la prima iterazione che vai a fare, dopo di che questo elemento, come detto anche da Wikipedia, andra' ad annullarsi.

Per prova, fai girare il tuo algoritmo con qualche scritta di debug:


void printMat (float a[4][4])
{
printf ("Input matrix\n");
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
printf ("%f ", a[i][j]);
}
printf ("\n");
}
}

void inversa(float a[4][4])
{
printMat (a);
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++)
{
printf ("PIVOT (%d,%d): calculated from -((%f)/(%f) {a[%d][%d] / a[%d][%d]) = ", z, i, a[i+1][z], a[z][z], i+1, z, z, z);
c = -(a[i+1][z] / a[z][z]);
printf ("%f\n", c);
for (j=z;j<4;j++)
{
printf (" a[i][j] = c * a[i][j]; {i=%d j=%d c=%f a[i][j]=%f} RESULT: ", i, j, c, a[i][j]);
a[i][j]= c * a[i][j];
printf ("%f\n", a[i][j]);

printf (" a[i + 1][j] = a[i][j] + a[i + 1][j] {a[i][j]=%f a[i + 1][j]=%f} RESULT: ", a[i][j], a[i+1][j]);
a[i + 1][j] = a[i][j] + a[i + 1][j];
printf ("%f\n", a[i+1][j]);

printf (" b[i][j]= c * b[i][j] {c=%f b[i][j]=%f} RESULT: ", c, b[i][j]);
b[i][j]= c * b[i][j];
printf ("%f\n", b[i][j]);

printf (" b[i + 1][j] =b[i][j]+ b[i + 1][j]; {b[i][j]=%f b[i + 1][j]=%f} RESULT:", b[i][j], b[i+1][j]);
b[i + 1][j] =b[i][j]+ b[i + 1][j];
printf ("%f\n", b[i+1][j]);
}
printf ("After cycle, a[z=%d+1][z=%d+1] = %f\n", z, z, a[z+1][z+1]);
}
}

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];

}}}

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]);
}
}
}


(ovviamente il primo blocco, quello con tutte le scritte di debug, e' importante per questa investigazione).

Vedrai che l'elemento e' andato a zero e, se non ricordo male, l'algoritmo e' stato fatto per funzionare cosi'.

Penso che sia necessario scegliere il pivot fra gli elementi non nulli, scandendoli sulla riga prescelta (come d'altronde hai riportato da Wikipedia)

[Thunderx]

Hai ragione, ma probabilmente non ci siamo capiti.
Quando devi calcolare c, utilizzi un elemento sulla diagonale per fare la divisione, confidente del fatto che la tua matrice di test ha la diagonale non nulla:

c = -(a[i+1][z] / a[z][z]);

In realta', questo e' vero per la prima iterazione che vai a fare, dopo di che questo elemento, come detto anche da Wikipedia, andra' ad annullarsi.

Per prova, fai girare il tuo algoritmo con qualche scritta di debug:


void printMat (float a[4][4])
{
printf ("Input matrix\n");
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
printf ("%f ", a[i][j]);
}
printf ("\n");
}
}

void inversa(float a[4][4])
{
printMat (a);
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++)
{
printf ("PIVOT (%d,%d): calculated from -((%f)/(%f) {a[%d][%d] / a[%d][%d]) = ", z, i, a[i+1][z], a[z][z], i+1, z, z, z);
c = -(a[i+1][z] / a[z][z]);
printf ("%f\n", c);
for (j=z;j<4;j++)
{
printf (" a[i][j] = c * a[i][j]; {i=%d j=%d c=%f a[i][j]=%f} RESULT: ", i, j, c, a[i][j]);
a[i][j]= c * a[i][j];
printf ("%f\n", a[i][j]);

printf (" a[i + 1][j] = a[i][j] + a[i + 1][j] {a[i][j]=%f a[i + 1][j]=%f} RESULT: ", a[i][j], a[i+1][j]);
a[i + 1][j] = a[i][j] + a[i + 1][j];
printf ("%f\n", a[i+1][j]);

printf (" b[i][j]= c * b[i][j] {c=%f b[i][j]=%f} RESULT: ", c, b[i][j]);
b[i][j]= c * b[i][j];
printf ("%f\n", b[i][j]);

printf (" b[i + 1][j] =b[i][j]+ b[i + 1][j]; {b[i][j]=%f b[i + 1][j]=%f} RESULT:", b[i][j], b[i+1][j]);
b[i + 1][j] =b[i][j]+ b[i + 1][j];
printf ("%f\n", b[i+1][j]);
}
printf ("After cycle, a[z=%d+1][z=%d+1] = %f\n", z, z, a[z+1][z+1]);
}
}

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];

}}}

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]);
}
}
}


(ovviamente il primo blocco, quello con tutte le scritte di debug, e' importante per questa investigazione).

Vedrai che l'elemento e' andato a zero e, se non ricordo male, l'algoritmo e' stato fatto per funzionare cosi'.

Penso che sia necessario scegliere il pivot fra gli elementi non nulli, scandendoli sulla riga prescelta (come d'altronde hai riportato da Wikipedia)
si hai ragione in quanto il primo coefficieNTE moltiplica la prima riga per 0 ed ecco il patatrac.quindi dovrei fare in modo di non far eseguire iil ciclo se a[i+1][j] è uguale a 0
che ne dici di questa soluzione?
for (z=0;z<3;z++) {
for (i=z;i<3;i++) {
if (a[i+1][j]=0){ break;}
else {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];
}
}
}

[sottovento]

si hai ragione in quanto il primo coefficieNTE moltiplica la prima riga per 0 ed ecco il patatrac.quindi dovrei fare in modo di non far eseguire iil ciclo se a[i+1][j] è uguale a 0
che ne dici di questa soluzione?
for (z=0;z<3;z++) {
for (i=z;i<3;i++) {
if (a[i+1][j]=0){ break;}
else {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];
}
}
}

Ciao
purtroppo non sono un esperto di questo algoritmo.
Quello che posso dirti e' che c'e' una svista: hai scritto

if (a[i+1][j]=0){ break;}

(che assegna 0 all'elemento a[i+1][j])
Immagino volessi scrivere:

if (a[i+1][j]==0){ break;}


A dirla tutta non la capisco bene, poiche' la tua divisione era

c = -(a[i+1][z] / a[z][z]);

per cui mi sarei aspettato un controllo su a[z][z].

Comunque, leggendo wikipedia, sembra che non debba saltare il ciclo nel caso l'elemento sia zero, ma che si debba semplicemente scegliere un altro elemento come pivot (ovviamente se esiste). Nel caso un elemento diverso da zero non esista, allora la riga in oggetto deve essere eliminata (in tal caso otterrai un rango inferiore). Se comunque parti da una matrice per la quale hai gia' verificato che esiste l'inversa, questo non dovrebbe capitarti.

Scusa la risposta stringata ma sono ancora al lavoro.

[Thunderx]

Ciao
purtroppo non sono un esperto di questo algoritmo.
Quello che posso dirti e' che c'e' una svista: hai scritto

if (a[i+1][j]=0){ break;}

(che assegna 0 all'elemento a[i+1][j])
Immagino volessi scrivere:

if (a[i+1][j]==0){ break;}


A dirla tutta non la capisco bene, poiche' la tua divisione era

c = -(a[i+1][z] / a[z][z]);

per cui mi sarei aspettato un controllo su a[z][z].

Comunque, leggendo wikipedia, sembra che non debba saltare il ciclo nel caso l'elemento sia zero, ma che si debba semplicemente scegliere un altro elemento come pivot (ovviamente se esiste). Nel caso un elemento diverso da zero non esista, allora la riga in oggetto deve essere eliminata (in tal caso otterrai un rango inferiore). Se comunque parti da una matrice per la quale hai gia' verificato che esiste l'inversa, questo non dovrebbe capitarti.

Scusa la risposta stringata ma sono ancora al lavoro.

Ciao sotto vento....ho individuato l'errore .....era un indice messo male nell'algoritmo (d)......alla fine tutto apposto, ti ringrazio per l'aiuto.
Un saluto dall'italia