[C++] Programma di calcolo matrice inversa

[xaler]

Salve a tutti
Stavo scrivendo un programma per calcolare la matrice inversa di una matrice quadrata di qualsiasi ordine e avevo un paio di dubbi/problemi.
Spiego brevemente il funzionamento generale del programma:
Viene chiesto all'utente di inserire l'ordine della matrice e dunque viene creata una matrice di tale ordine (M) che viene poi inizializzata tramite una serie di cout e cin.
Viene poi creata un'altra matrice (m) con lo stesso numero di righe ma doppio numero di colonne. Nella prima metà delle colonne verrà inserita la matrice creata in precedenza (M) mentre nella seconda metà viene inserità la matrice identità.
Quindi tramite eliminazione di gauss (sopra e sotto la diagonale) e tramite la divisione per i pivot ottengo la mia matrice inversa. Ci sono un paio di problemi però:
1) prima di tutto come faccio a deallocare completamente tutto lo spazio in memoria occupato ? devo distruggere ad uno ad uno tutti i puntatori dinamici o c'è un modo per liberare la memoria occupata direttamente ?
2) quando faccio l'eliminazione di gauss sotto la diagonale (vedere sui commenti al codice) metto l'if(m[i][i]=0) per dire che se c'è un pivot che è 0 voglio che questa riga si sposti verso il basso e la riga sottostante salga verso sopra in modo da avere sempre in diagonale elementi non nulli. L'errore riscontrato è una serie di nan e inf . La parte di codice incriminato è

double *tmp; 
tmp=new double[2*o];
for(int j=0;j<o-1;j++)
{
      for(int i=j+1;i<o;i++)
      {  
      if(m[i][j]!=0) 
      {
      double mol=m[i][j]/m[j][j];
      for(int k=0;k<2*o;k++)
      {
      tmp[k]=mol*m[j][k];
      };
      for(int k=0;k<2*o;k++) 
      {
      m[i][k]-=tmp[k];
      };
      }
      if(m[i][i]=0)
      {
      double *tmp2;
      double *tmp3;
      tmp2= new double[2*o];
      tmp3= new double[2*o];
      for(int k=0; k<2*o; k++) {tmp2[k]=m[i][k]; tmp3[k]=m[i+1][k];}
      for(int k=0; k<2*o; k++) {m[i][k]=tmp3[k]; m[i+1][k]=tmp2[k];}
      }
      }
}

Spero di non essere stato eccessivo nella richiesta !
Grazie anticipatamente
PS se serve tutto il codice lo posso postare .. è che sono 150 righe quindi non so quanto possa tornare comodo. Sono nelle vostre mani

[MegaAlchimista]

if(m[i][i]==0)

stai attento. Se metti if(m[i][i]=0) dai agli elementi il valore 0

[MItaly]

1) Dato che hai allocato una "finta matrice multidimensionale" in stile C, devi deallocare le righe una a una (e poi deallocare anche il vettore con i puntatori alle righe).

Nota, tra parentesi, che nei calcoli in virgola mobile non è quasi mai corretto confrontare direttamente con 0: dato che hanno precisione finita e nei conti si possono verificare errori di approssimazione quello che si fa in genere è verificare se sono sotto un certo epsilon.

[xaler]

Continua a non funzionare con "TUTTE" le matrici quadrate
se utilizzo questa matrice:
-2 2 -3
-1 1 3
2 0 -1

Il programma va in tilt e si riempie di nan e inf...
il mio codice è scritto sotto... è 3 giorni che lo scrivo (sono un neofita) e confido molto in voi

#include<iostream>
#include<math.h>
#include<cstdlib>
using namespace std;

// Alloco la memoria per una matrice di ordine "o" e la inizializzo
int funz1()
{   
int o;
cout << "Dammi l'ordine della tua matrice: ";
cin >> o;
double **M;
M= new double*[o];
for(int i = 0; i < o; i++)
   {
    M[i] = new double[o];       
    };
double val;
for(int k=0; k<o ; k++)
    {
       for(int l=0; l<o; l++)
       {
       cout << "Riga " << k+1 << " Colonna " << l+1 << ": " ;
       cin >> val;
       M[k][l]=val;
       };
    };
cout << "\nLa tua matrice è:\n" << endl;
for(int i=0; i<o; i++)
{ 
     for(int j=0; j<o; j++) 
     {
      cout << M[i][j] << "	";
     } 
cout << "\n" << endl;
}

//Alloco la memoria per una nuova matrice con lo stesso numero di righe ma doppio numero di colonne ( |A||I| )

double **m;
m= new double*[o];
for(int i=0; i<2*o; i++)
{
m[i]= new double[2*o];
}

for(int i=0; i<o; i++)
{
 for (int k=0; k<o; k++)
  {
  m[i][k]=M[i][k];
  }
}

for(int i=0; i<o; i++)
{
     for(int k=o; k<(2*o); k++)
     {
     if(i==(k-o)){m[i][k]=1;}
     else{m[i][k]=0;}
     };
}
cout << "Affianco la matrice identità:\n" << endl;
for(int i=0; i<o; i++)
{
   for(int k=0; k<2*o; k++)
   {
   cout << m[i][k] << "	";
   }
cout << "\n" << endl;
}

//Cominciano le operazioni di calcolo per la matrice inversa

double *tmp; 
tmp=new double[2*o];
for(int j=0;j<o-1;j++)
{
      for(int i=j+1;i<o;i++)
      {  
      if(m[i][j]!=0) 
      {
      double mol=m[i][j]/m[j][j];
      for(int k=0;k<2*o;k++)
      {
      tmp[k]=mol*m[j][k];
      };
      for(int k=0;k<2*o;k++) 
      {
      m[i][k]-=tmp[k];
      };
      }
      if(m[i][i]==0)
      {
      double *tmp2;
      double *tmp3;
      tmp2= new double[2*o];
      tmp3= new double[2*o];
      for(int k=0; k<2*o; k++) {tmp2[k]=m[i][k]; tmp3[k]=m[i+1][k];}
      for(int k=0; k<2*o; k++) {m[i][k]=tmp3[k]; m[i+1][k]=tmp2[k];}
      }
      }
}

cout << "Faccio l'eliminazione di Gauss sotto la diagonale:\n" << endl;
for(int i=0; i<o; i++)
{
   for(int k=0; k<2*o; k++)
   {
   cout << m[i][k] << "	";
   };
cout << "\n" << endl;
};

for(int j=o-1;j>0;j--)
{  
    for(int i=j-1;i>=0;i--)
      {
         if(m[i][j]!=0) 
        {
        double mol=m[i][j]/m[j][j];
        for(int k=0;k<2*o;k++) {tmp[k]=mol*m[j][k];};
        for(int k=0;k<2*o;k++) {m[i][k]-=tmp[k];};
        }
      };
};

//Controllo del determinante
double d=m[0][0];
for(int i=0; i<o; i++)
{d=d*m[i][i];}
cout << "il determinante è: " << d << "\n" << endl;
if(d=0){cout << " il determinante è 0. Non è possibile calcolare l'inversa di questa matrice ! " << endl; exit(0);}

cout << "Faccio l'eliminazione di Gauss sopra la diagonale:\n " << endl;
for(int i=0; i<o; i++)
{
   for(int k=0; k<2*o; k++)
   {
   cout << m[i][k] << "	";
   };
cout << "\n" << endl;
};

for(int i=0;i<o;i++)
{
    if(m[i][i]!=1) 
    {
    double mol=m[i][i];
      for(int k=0;k<2*o;k++)
      {
      m[i][k]=m[i][k]/mol;
      };
    }
};

cout << "Divido per i pivots:\n" << endl;
for(int i=0; i<o; i++)
{
   for(int k=0; k<2*o; k++)
   {
   cout << m[i][k] << "	";
   };
cout << "\n" << endl;
};

} 

int main()
{
cout << "Benvenuto in MatrixInverter" << endl;
funz1();
}

Come potete vedere subito dopo il commento //cominciano le operazioni di calcolo della matrice inversa
li sto cercando di fare una eliminazione di gauss che mi permetta anche di scambiare la riga con quella sotto nel caso in cui un pivot diventi 0... una volta compilato ed eseguito si presentano 2 problemi eseguendo la matrice sopra scritta:
1) mi si riempie di nan e inf vari come se avesse utilizzato uno 0 per le proprie operazioni (e non capisco a quale 0 possa poggiarsi.
2) le righe vengono invertite ma ovviamente quella sotto che viene portata su "aggira" il for e non viene ridotta dalle mosse di gauss

come posso cavarci fuori qualcosa ? premetto che ci ho pensato molto ma non giungo a conclusioni esaudienti

[MegaAlchimista]

Purtroppo ora non posso provare il codice quindi non posso aiutarti più di tanto, comunque dove calcoli il determinante hai di nuovo scritto

if(d=0)
cout << "il determinante è 0";

per ora posso solo dirti di ricontrollare questa cosa, e fare come suggerito da MItaly e non confrontare i double/float con 0 , ma per esempio assicurati che il loro valore assoluto sia minore di 0.000001.

Se dopo questi riadattamenti non va, metti quella matrice che ti crea problemi ed esegui il codice in debug, e segui i risultati di tutti i calcoli finchê non trovi la riga nella quale si insidia l'errore

[xaler]

Il problema è il seguente:
quando nella diagonale compaiono degli 0 (prima o dopo l'eliminazione di Gauss) questi vengono usati per ridurre la matrice a forma triangolare e ovviamente viene diviso tutto per 0 e il programma restituisce nan o inf. Allora ho creato la condizione con l'if che se si trova un elemento in diagonale con 0 (quelli con indici uguali dunque) l'intera riga viene spostata verso il basso e quella sotto viene spostata su. Il problema è che la riga che da sotto passa sopra bypassa l'eliminazione di gauss e non viene ridotta. Avete idee ?


In piu
come faccio a deallocare tutto ?
alla fine della funzione io scrivo:

for(int i=0;i<o;i++) {delete M[i];delete m[i];};
for(int i=o;i<2*o;i++) {delete m[i];};
delete m;
delete M;
delete tmp;

Ma mi viene restituito un errore in fase di esecuzione:

*** glibc detected *** ./a.out: double free or corruption (out): 0x0000000000be20e0 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x7e626)[0x7f5b38d98626]
./a.out[0x401a80]
./a.out[0x401b2e]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed)[0x7f5b38d3b76d]
./a.out[0x400939]
======= Memory map: ========
00400000-00402000 r-xp 00000000 08:05 786731                             /home/dario/a.out
00602000-00603000 r--p 00002000 08:05 786731                             /home/dario/a.out
00603000-00604000 rw-p 00003000 08:05 786731                             /home/dario/a.out
00be2000-00c03000 rw-p 00000000 00:00 0                                  [heap]
7f5b38806000-7f5b3880a000 rw-p 00000000 00:00 0 
7f5b3880a000-7f5b3881f000 r-xp 00000000 08:05 1445543                    /lib/x86_64-linux-gnu/libgcc_s.so.1
7f5b3881f000-7f5b38a1e000 ---p 00015000 08:05 1445543                    /lib/x86_64-linux-gnu/libgcc_s.so.1
7f5b38a1e000-7f5b38a1f000 r--p 00014000 08:05 1445543                    /lib/x86_64-linux-gnu/libgcc_s.so.1
7f5b38a1f000-7f5b38a20000 rw-p 00015000 08:05 1445543                    /lib/x86_64-linux-gnu/libgcc_s.so.1
7f5b38a20000-7f5b38b19000 r-xp 00000000 08:05 1445554                    /lib/x86_64-linux-gnu/libm-2.15.so
7f5b38b19000-7f5b38d18000 ---p 000f9000 08:05 1445554                    /lib/x86_64-linux-gnu/libm-2.15.so
7f5b38d18000-7f5b38d19000 r--p 000f8000 08:05 1445554                    /lib/x86_64-linux-gnu/libm-2.15.so
7f5b38d19000-7f5b38d1a000 rw-p 000f9000 08:05 1445554                    /lib/x86_64-linux-gnu/libm-2.15.so
7f5b38d1a000-7f5b38ecd000 r-xp 00000000 08:05 1445522                    /lib/x86_64-linux-gnu/libc-2.15.so
7f5b38ecd000-7f5b390cc000 ---p 001b3000 08:05 1445522                    /lib/x86_64-linux-gnu/libc-2.15.so
7f5b390cc000-7f5b390d0000 r--p 001b2000 08:05 1445522                    /lib/x86_64-linux-gnu/libc-2.15.so
7f5b390d0000-7f5b390d2000 rw-p 001b6000 08:05 1445522                    /lib/x86_64-linux-gnu/libc-2.15.so
7f5b390d2000-7f5b390d7000 rw-p 00000000 00:00 0 
7f5b390d7000-7f5b391b9000 r-xp 00000000 08:05 926486                     /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f5b391b9000-7f5b393b8000 ---p 000e2000 08:05 926486                     /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f5b393b8000-7f5b393c0000 r--p 000e1000 08:05 926486                     /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f5b393c0000-7f5b393c2000 rw-p 000e9000 08:05 926486                     /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f5b393c2000-7f5b393d7000 rw-p 00000000 00:00 0 
7f5b393d7000-7f5b393f9000 r-xp 00000000 08:05 1445502                    /lib/x86_64-linux-gnu/ld-2.15.so
7f5b395db000-7f5b395dc000 rw-p 00000000 00:00 0 
7f5b395f4000-7f5b395f9000 rw-p 00000000 00:00 0 
7f5b395f9000-7f5b395fa000 r--p 00022000 08:05 1445502                    /lib/x86_64-linux-gnu/ld-2.15.so
7f5b395fa000-7f5b395fc000 rw-p 00023000 08:05 1445502                    /lib/x86_64-linux-gnu/ld-2.15.so
7fff67cf2000-7fff67d13000 rw-p 00000000 00:00 0                          [stack]
7fff67dff000-7fff67e00000 r-xp 00000000 00:00 0                          [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]
Aborted

Si può risolvere ?

[xaler]

UP !
sto impazzendo !

[xaler]

UP !

[MItaly]

Trattandosi di array, non di deallocano con `delete`, ma con `delete []`:

for(int i=0;i<o;i++) {delete[] M[i];delete[] m[i];};
for(int i=o;i<2*o;i++) {delete[] m[i];};
delete[] m;
delete[] M;
delete[] tmp;

(anche se alcuni compilatori tendono a perdonare questo errore)