[C++] Potenza di Matrice

[Mancio1]

Ciao a tutti, chiedo aiuto in questo forum per risolvere un esercizio proposto nelle dispense dal mio professore di informatica. É un esercizio riguardante le potenze di matrici in c++ (come da titolo). Premetto che, seppur ho passato l'appello di Geometria ed Algebra Lineare, non mi sono mai e poi mai trovato nella situazione di dover calcolare una potenza di matrice, pertanto sono nel buio totale. L'unica cosa che mi viene in mente, sarebbe una moltiplicazione per se stessa della matrice, per un determinato numero di volte...quindi pensavo ad un algoritmo che mi permettesse questo. Vi posto intanto l'intestazione del problema, ed il codice fin'ora scritto:

Scrivere un programma C++ che legge da tastiera (e da file) una matrice m di dimensioni N*M
calcola la matrice trasposta di m, mT
calcola la matrice prodotto m* mT
Per una buona organizzazione del codice vanno definiti e usati i seguenti sottoprogrammi
void calcolaTrasposta(int m[N][M], int t[M][N])
void calcolaProdotto(int m1[N][M], int m2[M][N], int prodotto[N][N])
__________________________________________________________________________________________________
Aggiungere al programma precedente le funzionalità necessarie per calcolare l'n-esima potenza del prodotto di un matrice m per la sua trasposta, ovvero ( m*mT )n
Per una buona organizzazione del codice vanno aggiunti, ed usati almeno i seguenti sottoprogrammi
void calcolaPotenza(int m[N][N], int n, int potenza[N][N])
che calcola mn e restituisce il risultato in potenza
void unita(int m[N][N])
che inizializza m alla matrice unità di dimensione N
void copiaMatrice(int s[N][N], int d[N][N])
che copia elemento per elemento la matrice s nella matrice d

P.S. sinceramente questa volta, non riesco a capire nemmeno perché mi abbia dato questi ultimi sottoprogrammi...

_________________________________________________________________________________________________
Questo é il codice, che riguarda solo la prima parte (perfettamente funzionante )

#include <iostream>
using namespace std;

int const N=3;
int const M=4;

void leggiMatrice(int matrice[N][M]);
void stampaMatrice(int matrice[N][M]);
void calcolaTrasposta(int m[N][M], int t[M][N]);
void calcolaProdotto(int m1[N][M],int m2[M][N], int prodotto[N][N]);
void calcolaPotenza(int m[N][N], int n, int potenza[N][N]);


int main(){
    int matrice[N][M], trasposta[M][N],prodotto[N][N];
    leggiMatrice(matrice);
    stampaMatrice(matrice);
    calcolaTrasposta(matrice,trasposta);
    calcolaProdotto(matrice,trasposta,prodotto);
system("PAUSE");
}

void leggiMatrice(int matrice[N][M]){
    cout << "Inserisci la matrice " << N << " X " << M << ": " <<endl;
    for(int i=0;i<N;i++){
        for(int j=0;j<M;j++){
            cin >> matrice[i][j];
        }
    }
    return;
            
}

void stampaMatrice(int matrice[N][M]){
    cout << "La matrice immessa e': " <<endl;
    for(int i=0;i<N;i++){
        for(int j=0;j<M;j++){
            cout << matrice[i][j] << "\t";
        }
    cout <<endl;
    }
    return;

}
void calcolaTrasposta(int m[N][M], int t[M][N]){
    cout << "La matrice trasposta e': " <<endl;
    for(int i=0;i<N;i++){
        for(int j=0;j<M;j++){
            t[j][i]=m[i][j]; //scopri perché solo in questo ordine funziona, e non scrivendo "m[i][j]=t[j][i]"
        }
    }
    for(int i=0;i<M;i++){
         for(int j=0;j<N;j++){
             cout << t[i][j] << "\t";
         }
         cout <<endl;
     }
     return;    
}
void calcolaProdotto(int m1[N][M],int m2[M][N], int prodotto[N][N]){
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            prodotto[i][j] = 0; // inizializziamo a zero tutti gli elementi della matrice, grazie ai due cicli "for" nidificati
            for(int k=0;k<M;k++){ //(spiegazione accurata) qui usiamo un parametro qualsiasi per bloccare una volta la "riga", una volta una "colonna"
                prodotto[i][j] += m1[i][k]*m2[k][j]; // usiamo "+=" per andare a sostituire lo "0" con il prodotto riga per colonna
            }
        }
    }
    cout << "La matrice prodotto delle due e': " <<endl;
    for (int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            cout << prodotto[i][j] << "\t";
        }
        cout <<endl;
    }
    return;
}

void calcolaPotenza(int m[N][N], int n, int potenza[N][N]){
    
}

[apatriarca]

Ma no.. Implementando la potenza come un ciclo di prodotti come hai fatto prima.. Però stavolta invece di fare il prodotto richiami la funzione..

[vict85]

In quella struttura non c'é nulla di “ben organizzato”, almeno a mio avviso. Sta chiedendo di fare la moltiplicazione di una matrice per la sua trasposta. Ora, non per fare il pignolo, ma calcolare la trasposta in questo caso è semplicemente stupido: è uno spreco di memoria e tempo computazionale notevole. Mi spiego meglio, il normale prodotto matriciale è una moltiplicazione riga per colonna, il prodotto di una matrice per una trasposta è un prodotto riga per riga (della matrice non trasposta). Nel caso in cui stai facendo questa operazione sulla stessa matrice allora tu fai inutilmente una copia della matrice di partenza e ne calcoli inutilmente la trasposta. Senza ridurre minimamente la complessità della moltiplicazione, anzi, finendo addirittura per mancare di sfruttare il fatto che la moltiplicazione per una trasposta sia notevolmente cache friendly (mentre la normale moltiplicazione con i cicli strutturati in maniera classica è l'apoteosi dei cache missing). Insomma penso che sia una struttura notevolmente più lenta (penso nell'ordine di 3-4 volte almeno sull'ottimale a singolo processo).

Detto questo, nel mio primo corso di informatica un algoritmo del genere dell'esponenziazione veloce non era considerato fuori-programma. Anche perché è un esempio piuttosto classico di funzione ricorsiva, e i vantaggi rispetto al metodo naive sono tanti, specialmente quando la moltiplicazione ha una complessità cubica. Penso che provare a implementarlo possa essere educativo.

Per concludere la mia disquisizione sulle matrici, ti faccio notare che nel ciclo:

for(i=0; i<m; ++i)
  for(j=0; j<n; ++j)
    for(k=0; k<s; ++k)
      c[i][j] += a[i][k] * b[k][j];

possiede la interessante proprietà che gli indici possono ‘permutare’, cioé esistono varianti in cui l'ordine di ijk è differente. Le performance dei vari cicli possono essere notevolmente diverse anche se ad un programmatore poco esperto può apparire poco intuitivo.
Nel tuo caso particolare, se non ho sbagliato i calcoli, il migliore, a meno di metodi più complessi, dovrebbe essere il ciclo ikj cioé il seguente:

for(i=0; i<m; ++i)
  for(k=0; k<s; ++k)
    for(j=0; j<n; ++j)
      c[i][j] += a[i][k] * b[k][j];

Affinché quei cicli funzionino devi aver azzerato la matrice c. In caso contrario devi adattare adeguatamente quei cicli.

[Mancio1]

Ho letto quello che mi hai scritto, ho afferrato il concetto della velocitá dei cicli, e se dovessimo parlare di "ben organizzato" ti direi, lasciamo perdere .

Quello che non mi é chiaro é come usare la funzione nel ciclo di prodotti. Cioé, se la funzione é una void, e quindi una procedura, non mi ritornerá alcun valore. Come la inserisco in calcolaPotenza?
Scusa se ti faccio domande banali, peró non credo di aver mai affrontato questi tipo di problem all'universitá durante il corso, diciamo che era tutto piú facile.

[vict85]

Non penso tu abbia compreso sinceramente, anche perché erano considerazioni un po' tecniche (per capire perché un certo ciclo è più veloce di un altro devi aver presente come la cpu lavora sulla memoria e che operazioni sono più veloci). Inoltre se quella del professore era una struttura criticabile per il suo obiettivo finale, la tua è carente dove invece la sua è buona. La sua struttura infatti ha una funzione per ogni singola operazione, la tua invece dimostra ancora un modo di vedere le funzioni decisamente da principiante (la colpa è in gran parte di libri e professori che spiegano le cose male).

In particolare metti dei cout e dei cin dove non sono necessari. Questo è un gravissimo errore di design che devi perdere in fretta se vuoi occuparti di programmazione in futuro.

Prima però di trattare aspetti di design ti spiego qualcosa che tu non hai capito: il fatto che una funzione non ritorni nulla non significa affatto che non modifichi nulla. In particolare, in questo caso, modifica le matrici che tu passi alla funzione. Il fatto è che quando passi un array, in C, lo passi sempre per riferimento o in altre parole è come se tu avessi un puntatore. In pratica quindi quando lavori all'interno di quelle funzioni modifichi le matrici che hai definito nel main.
Per quanto riguarda quindi calcola potenza il codice che devi usare è qualcosa di questo tipo:

unita(prodotto);
for(int i = 0; i<n; ++i)
  calcolaProdotto(prodotto, m, prodotto);

In particolare NON devi richiedere n nella funzione, né visualizzare alcunché a video. Questa è comunque la versione naif. Ci sono inoltre altri modi in cui una funzione che non ritorna nulla possa modificare qualcosa al di fuori di se stessa, ma per i tuoi scopi puoi ignorarli.

Riguardo al design devi toglierti la mania di visualizzare a video che cosa fanno le funzioni. Una funzione deve fare una sola cosa e nient'altro. Se deve calcolare il prodotto allora deve calcolare un prodotto e non dire che lo deve fare e chiedere che elementi deve moltiplicare. Ci saranno altre parti del programma che se ne preoccuperanno anche perché spesso una funzione viene usata da altre funzioni e non direttamente dall'utente. O ancora non tutti i programmi usano la console. Senza considerare che se hai bisogno di un calcolo e ne devi fare 1000 intermedi tu non visualizzi tutti e mille i calcoli intermedi. In pratica se vuoi mettere quei cout/cin, li devi mettere nel main.

[Mancio1]

Ragazzi vi ringrazio delle risposte, cercheró di apprendere al meglio I vostri consigli, perché sí, seppur non fa parte del corso di Ing. Meccanica, mi piace molto come studio, e soprattutto, credo sia molto utile se non indispensabile ai giorni nostri. Ci sto provando da ieri, ancora non riesco a tirare fuori questo benedetto codice, ho l'appello a breve (anche se non capiterá mai nulla del genere) vi chiedo gentilmente di postarmelo fatto cosicché lo possa confrontare, studiare e capirne gli errori commessi. Grazie

[vict85]

Provandolo a fare ho trovato un errore, ritengo piuttosto grave e problematico, del design del tuo professore. Se la tua funzione calcolaPotenza è scritta per matrici m1[N][M] e m2[M][N] come posso usarla per matrici m1[N][N] e m2[N][N]? La risposta è che il compilatore ti da errore se ci provi. La soluzione più semplice è scrivere una seconda funzione o fare come hai fatto tu, ma entrambe le opzioni sono poco pratiche e cattivi design. Il problema è che di metodi elementari per risolvere il problema non ne conosco. Il più semplice, nonché quello usato abitualmente è quello che ho usato io:

#include <iostream>
using namespace std;


// utility definite da te
void leggiMatrice(int * const matrice, unsigned const Idim, unsigned const Jdim);
void stampaMatrice(int const * const matrice, unsigned const Idim, unsigned const Jdim);

// funzioni definite dal professore
void calcolaTrasposta(int const * const m, int * const t, unsigned const Idim, unsigned const Jdim);
void calcolaProdotto(int const * const m1, int const * const m2, int * const prodotto, unsigned const Idim, unsigned const Jdim, unsigned const Kdim);
void calcolaPotenza(int const * const m, int * const potenza, unsigned const esponente, unsigned const N);
void unita(int * const m, unsigned const N);
void copiaMatrice(int const * const s, int * const d, unsigned const Idim, unsigned const Jdim);


int main()
{
    int const N=3;
    int const M=4;
    int matrice[N][M], trasposta[M][N], prodotto[N][N], n, potenza[N][N];
    
    cout << "Inserisci la matrice " << N << " x " << M << ": " <<endl;
    leggiMatrice(&matrice[0][0], N, M);
    cout << endl;
    cout << "La matrice da te inserita e' :" << endl;
    stampaMatrice(&matrice[0][0], N, M);
    
    cout << "Sto calcolando la trasposta... ";
    calcolaTrasposta(&matrice[0][0], &trasposta[0][0], N, M);
    cout << "Il risultato ottenuto e' :" << endl;
    stampaMatrice(&trasposta[0][0], M, N);
    
    cout << "Sto calcolando il prodotto tra la matrice immessa e la sua trasposta... ";
    calcolaProdotto(&matrice[0][0],&trasposta[0][0],&prodotto[0][0], N, N, M);
    cout << "Il risultato ottenuto e' :" << endl;
    stampaMatrice(&prodotto[0][0], N, N);
    
    cout << "A che valore vuoi elevare potenza? ";
    cin >> n;
    cout << "Sto elevando il prodotto alla potenza " << n << "-esima... ";
    calcolaPotenza(&prodotto[0][0], &potenza[0][0], n, N);
    cout << "Il risultato ottenuto e' :" << endl;
    stampaMatrice(&potenza[0][0], N, N);
    
    //system("PAUSE"); //questo non è necessario se fai partire il programma da console o usando ctrl+F5
}

void leggiMatrice(int * const matrice, 
                  unsigned const Idim, unsigned const Jdim)
{
    
    for(unsigned i=0; i != Idim; ++i) {
        for(unsigned j=0; j != Jdim; ++j) {
            cin >> matrice[i*Jdim + j];
        }
    }
}

void stampaMatrice(int const * const matrice,
                   unsigned const Idim, unsigned const Jdim)
{
    for(unsigned i=0; i != Idim; ++i) {
        for(unsigned j=0; j != Jdim; ++j) {
            cout << matrice[i*Jdim + j] << "\t";
        }
        cout <<endl;
    }
}

void calcolaTrasposta(int const * const m, 
                      int       * const t,
                      unsigned const Idim, unsigned const Jdim)
{
    for(unsigned i=0; i != Idim; ++i)
        for(unsigned j=0; j != Jdim; ++j)
            t[j*Idim + i] = m[i*Jdim + j];
}


/*
 *  Attenzione questa funzione richiede che prodotto NON sia uguale a m1 o m2. In caso contrario è necessario fare prima una copia.
 *  Per risolvere il problema sarebbe necessario una analisi più ragionata dell'algoritmo. Ho usato un ciclo ikj.
 */
void calcolaProdotto(int const * const m1,
                     int const * const m2,
                     int       * const prodotto,
                     unsigned const Idim, unsigned const Jdim, unsigned const Kdim)
{
    for(unsigned i=0; i != Idim; ++i) {
        for(unsigned j=0; j != Jdim; ++j)
            prodotto[i*Jdim + j] = m1[i*Kdim] * m2[j];
        for(unsigned k = 1; k != Kdim; ++k)
            for(unsigned j=0; j != Jdim; ++j) 
                prodotto[i*Jdim + j] += m1[i*Kdim + k] * m2[k*Jdim + j];         
    }
}

/*
 *  Attenzione questa funzione richiede che potenza NON sia uguale a m. In caso contrario è necessario fare prima una copia.
 *  Per risolvere il problema sarebbe necessario una analisi più ragionata dell'algoritmo.
 */
void calcolaPotenza(int const * const m,
                    int       * const potenza,
                    unsigned const esponente, unsigned const N)
{
    unita(potenza, N);
    int * temp = new int[N*N];
    for(unsigned i=0; i != esponente; ++i) {
        copiaMatrice(potenza, temp, N, N); // la copia è necessaria perché il prodotto di matrici altrimenti non funziona.
        calcolaProdotto(temp, m, potenza, N, N, N);
    }
    delete[] temp;
}

void unita(int * const m, 
           unsigned const N)
{
    for(unsigned i = 0; i != N; ++i)
        for(unsigned j = 0; j != N; ++j)
            m[i*N + j] = (i==j)? 1 : 0;
}

void copiaMatrice(int const * const s,
                  int       * const d,
                  unsigned const Idim, unsigned const Jdim)
{
    for(unsigned i=0; i != Idim; ++i)
        for(unsigned j=0; j != Jdim; ++j)
            d[i*Jdim + j] = s[i*Jdim + j];
}

Di fatto quelle funzioni possono essere utilizzate anche per matrici dinamiche (anzi di fatto sono pensate con quello in mente).

Mi sono comunque accorto di aver fatto un errore prima. Nella potenza e nel prodotto è indispensabile che le varie matrici non siano la stessa matrice (altrimenti sovrascrivi la stessa locazione di memoria). Quindi nel calcolo della potenza ho dovuto creare una variabile temporanea per gestire questo problema. Nel design del professore questo aspetto era stato preso in considerazione (peccato che poi non ha fatto caso alle dimensioni).

[Mancio1]

Ringrazio per la pazienza e le spiegazioni. Devo studiare in maniera piú approfondita la materia, perché noto di non essere in grado di leggere il codice che mi hai scritto, non conosco proprio i comandi che sono stati usati. Sono peró soddisfatto perché ho imparato sicuramente qualcosa postando qui la mia domanda, spero di riuscire a passare l'appello di questa sessione estiva, e di continuare poi da autodidatta, magari postando altri quesiti qui o seguendo qualche corso. A presto!

[vict85]

Che comandi non capisci? In realtà non ci sono cose difficili da spiegare singolarmente. Hai fatto i puntatori?

[Mancio1]

No...guarda, il programma del professore é questo:
https://www.sites.google.com/site/roma3 ... -didattico

nel tuo codice cio che non ho ancora mai incontrato sono ad esempio:

"Unsigned"," m[i*N + j] = (i==j)? 1 : 0;", "&" (prima di un array), "*" (nel sottoprogramma)....cose cosí

[vict85]

&var[j] restituisce l'indirizzo dell'elemento i-j dell'array var. Serve nel caso tu voglia fare un puntatore a quell'elemento. Ma se non sai cos'è un puntatore allora ha poco senso spiegartelo.

Per quanto invece riguarda

var = condizione ? A : B; 

questo equivale a

if(condizione)
  var = A;
else
  var = B;

Si chiama operatore ternario.

Sinceramente penso che la cosa migliore sia spiegare il problema al professore e chiedere che fare. Il metodo probabilmente più semplice consiste nell'usare matrici quadrate. Se hai bisogno di maggiori spiegazioni su qual'è esattamente il problema te lo scrivo per bene.

[Mancio1]

Comodo quell'operatore ternario, non lo conoscevo, proveró ad usarlo piú spesso. Oggi ho il colloquio con il professore. Si scrivimi,che puó farmi molto comodo, grazie .

[vict85]

Il problema è che la funzione esplicita la dimensione degli array multidimensionali. Il C usa questa informazione per implementare correttamente l'operazione a[j]. Se quindi hai scritto una funzione per usare array multidimensionali di una certa dimensione non puoi usarlo per dimensioni diverse. Ad essere precisi devi aver esplicitato per bene ogni dimensione a parte quella più a destra. Ci sono delle precise ragioni legate alla memoria.

[Mancio1]

Perfetto, ho capito. Detto tra noi, se volessi che l'utente inserisse egli stesso la dimensione dell'array, dovrei usare il tipo unsigned che hai usato tu?

[vict85]

No, semplicemente ho usato unsigned perché i numeri in gioco erano non-negativi. Ma è una mia scelta stilistica.