#include <iostream>
using namespace std;
int main(){
// *** Allocazione della memoria *** //
double **A; int n;
cout << "\n\nCreazione matrice n x n\nInserire il valore di n: "; cin >> n; cout << "\n";
A=new double*[n];
for(int i=0;i<n;i++) A[i]=new double[n];
double **I;
I=new double*[n];
for(int i=0;i<n;i++) I[i]=new double[n];
// *** Popolazione della matrice identità *** //
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
if(i==j) I[i][j] = 1;
else I[i][j]=0;
// *** Inserimento degli elementi via tastiera *** //
cout << "\n************ INSERIMENTO *************\n";
for(int i=0;i<n;i++)
for(int j=0;j<n;j++){
cout << "Elem " << i << ", " << j << ": ";
cin >> A[i][j];
}
cout << "***************************************\n\n\n";
// *** Stampa della matrice *** //
cout << "\n***** MATRICE A ******\n\n";
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
cout << A[i][j] << "\t";
cout << "\n";
}
cout << "*********************\n\n";
// *** Costruzione della matrice [A|I] *** //
double **B;
B=new double*[n];
for(int i=0;i<n;i++) B[i]=new double[2*n];
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
B[i][j]=A[i][j];
int k=0;
for(int i=0;i<n;i++) {
for(int j=n;j<2*n;j++,k++)
B[i][j]=I[i][k];
k=0;
}
// *** Eliminazione sotto la diagonale principale *** //
double *tmp; tmp=new double[2*n];
for(int j=0;j<n-1;j++)
for(int i=j+1;i<n;i++)
if(B[i][j]!=0) {
double mol=B[i][j]/B[j][j];
for(int k=0;k<2*n;k++) tmp[k]=mol*B[j][k];
for(int k=0;k<2*n;k++) B[i][k]-=tmp[k];
}
// *** Eliminazione sopra la diagonale principale *** //
for(int j=n-1;j>0;j--)
for(int i=j-1;i>=0;i--)
if(B[i][j]!=0) {
double mol=B[i][j]/B[j][j];
for(int k=0;k<2*n;k++) tmp[k]=mol*B[j][k];
for(int k=0;k<2*n;k++) B[i][k]-=tmp[k];
}
// *** Ultimo step per ottenere la matrice a blocchi [I|A] *** //
for(int i=0;i<n;i++)
if(B[i][i]!=1) {
double mol=B[i][i];
for(int k=0;k<2*n;k++)
B[i][k]=B[i][k]/mol;
}
// *** Copia dell'inversa ottenuta *** //
double** Inv;
Inv=new double*[n];
for(int i=0;i<n;i++) Inv[i]=new double[n];
k=0;
for(int i=0;i<n;i++) {
for(int j=n;j<2*n;j++,k++)
Inv[i][k]=B[i][j];
k=0;
}
// *** Stampa della matrice Inversa*** //
cout << "\n***** MATRICE INVERSA ******\n\n";
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
cout << Inv[i][j] << "\t";
cout << "\n";
}
cout << "***************************\n\n";
// *** Azzeramento della matrice idendità utilizzata prima, che utilizzeremo in seguito come supporto *** //
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
I[i][j]=0;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(k=0;k<n;k++)
I[i][j]+=(A[i][k]*Inv[k][j]);
// *** Stampa della matrice identità ottenuta da A*A^-1 *** //
cout << "\n** Stampa di A*A^-1 **\n";
for(int i=0;i<n;i++) {
for(int j=0;j<n;j++)
cout << I[i][j] << "\t";
cout << "\n";
}
cout << "***********************\n";
// *** Deallocazione della memoria *** //
for(int i=0;i<n;i++) {delete I[i];delete A[i];}
delete I;
delete A;
for(int i=0;i<n;i++) delete B[i];
delete B;
delete tmp;
}
Tieni presente che per fare calcolo numerico c++ non è il massimo (per queste esigenze esistono altri software), dunque potresti spesso trovarti di fronte a numeri come 1.4322e-15 (ad esempio nell'ultima parte, il prodotto tra A e A^-1).
Chiaramente sono numeri minuscoli e capisci che lì ci dovrebbe essere uno zero.
Potresti anche fare un caso di test su questo "problema" per parlare della stabilità degli algoritmi, della cancellazione numerica ecc...
Il programma l'ho testato su linux e windows e sembra funzionare.
Tieni però presente che ci potrebbero volere dei controlli che io non ho fatto.
Due paroline sul codice:
- per calcolare A^-1 uso le operazioni elementare sulle righe della matrice, e sfrutto il fatto che data una matrice a blocchi
[A|I] mediante tali operazioni ottengo la matrice [I|A] (se hai bisogno, guardati qualcosa sull'eliminazione di gauss)
- L'allocazione avviene dinamicamente quindi alloco prima memoria per n puntatori e poi per ogni puntatore alloco n double
- Ho scritto tutto nel main per fare prima, chiaramente è da sistemare
- Non ho implementato la lettura da file, ma credo che la cosa più difficile era creare gli array e gestirli bene per calcolare l'inversa.
- non faccio controlli sui pivot durante l'eliminazione. Per ogni colonna dovresti cercare max(a_i) e scambiare le righe. Non l'ho fatto per sveltire il lavoro. Chiaramente se hai un pivot uguale a zero, il mio codice non funziona. Quindi ricorda di modificare quella parte.
- A deve essere non singolare, dunque se ti interessa fare quest'ulteriore controllo, devi calcolare il rango della matrice
- Ho fatto alcuni test con scilab e sembrano non esserci problemi, ovviamente fai più controlli che puoi perchè potrebbe essermi sfuggito qualcosa
Se non ti piace l'eliminazione, ricorda che c'è il metodo dei cofattori. Ti sconsiglio quello di risolvere il sistema n^2 naturalmente.
Non esistare a fare domande se hai dei dubbi