[C] matrice inversa 4x4 con algoritmo di Gauss

[sottovento]

Originariamente inviato da Thunderx

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]

Originariamente inviato da sottovento

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]

Originariamente inviato da Thunderx

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]

Originariamente inviato da sottovento

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