Newsgroup: it.comp.lang.c++
Titolo: Calcolo del determinante

Ho provato per esercizio a scrivere un programma che sia in grado di
calcolare il determinante di una matrice quadrata di qualsiasi ordine. Per
farlo ho utilizzato l'algoritmo di Laplace, cioč il determinante č pari alla
somma dei prodotti degli elementi di una riga per il complemento algebrico
dell'elemento stesso. Cosě ad esempio il det.della matrice
       | 1  3  2|
       | 0  2  3|  =  1* | 2  3 | - 3 * | 0  3 |  +  2 * | 0  2 | =  1 * ( 2
* | 1 | - 3 * | 1 |) -.....=
       | 2  1  1|           | 1  1 |          | 2  1 |            | 2  1 |
                     =  1 *(-1)  -  3 *(-6)  +  2 *(-4) = 9
Si tratta di un algoritmo ricorsivo in cui la condizione di stop si verifica
quando si ha una matrice del primo ordine.
Il problema č che, fissato come punto di riferimento il tempo per calcolare
una matrice 2 x 2, il calcolo di una 3x3 impiega circa 3 volte tanto, quello
di una 4x4 12 volte. In generale una n x n dovrebbe impiegare n! / 2 volte
una 2 x 2. Poichč una 10x10 sul mio PC impiega circa 35 secondi, una 20 x 20
dovrebbe restituire il risultato ai nostri discendenti visto che impiega piů
di 700.000 anni!
Non che mi serva attualmente fare calcoli di questo tipo, mi chiedevo
tuttavia se esistesse un altro algoritmo(magari non ricorsivo) che permetta
di velocizzare un po' i tempi.
Volendo invece mantenere questo algoritmo, cosa potrei cambiare nel mio
codice per aumentare le prestazioni? Io ho pensato che potrebbe essere utile
cercare la riga o la colonna che contiene il numero maggiore di zeri,
evitando in tal modo il calcolo, solo perň se ci sono degli zeri!
Spero di non avervi annoiato con queste mie curiositŕ.
Marco
#include <iostream>
#include <stdlib.h>
#include <iomanip>
#include <time.h>
typedef int Type;
class Matrix
{
friend ostream& operator<<(ostream&, const Matrix&);
friend int det(const Matrix&);
public:
     Matrix();
     Matrix(int, int);
     ~Matrix();
     Matrix(const Matrix&);
     void set(Type,int, int);
     void set_random();
     Matrix minore(int i,int j) const;
private:
     Type** m;
     int rows;
     int cols;
};
Matrix::Matrix() {
     m = 0;
     rows = 0;
     cols = 0;
}
Matrix::Matrix(int righe, int colonne) {
     rows = righe;
     cols = colonne;
     if(!(rows && cols)) {
          m = 0;
          return;
     }
     m = new Type*[rows];
     for(int i = 0; i < rows; i++)
          m[i] = new Type[cols];
     for(int i = 0; i < rows; i++)
          for(int j = 0; j < cols; j++)
               m[i][j] = 0;
}
Matrix::Matrix(const Matrix& M) {
     rows = M.rows;
     cols = M.cols;
     m = new Type*[rows];
     for(int i = 0; i < rows; i++)
          m[i] = new Type[cols];
     for(int i = 0; i < rows; i++)
          for(int j = 0;j < cols; j++)
               m[i][j] = M.m[i][j];
}
Matrix::~Matrix() {
     for(int i = 0; i < rows; i++)
          delete[] m[i];
     delete[] m;
}
void Matrix::set(Type elem, int i, int j) {
     if(i < 0 || i >= rows) {
          cout << "(" << i << "," << j << ")" << " out of range" << endl;
      return;
     }
     if(j < 0 || j >= cols) {
          cout << "(" << i << "," << j << ")" << " out of range" << endl;
      return;
     }
     m[i][j] = elem;
}
void Matrix::set_random() {
        srand(time(0));
     for(int i = 0; i < rows; i++)
          for(int j = 0; j < cols; j++)
               m[i][j] = rand()%4;
}
Matrix Matrix::minore(int noi, int noj) const{
     Matrix result(rows - 1, cols - 1);
     for(int si = 0, di = 0; si < rows; si++)
     {
          for(int sj = 0, dj = 0; sj < cols; sj++)
          {
           if((si != noi) && (sj != noj))
                result.set(m[si][sj], di, dj);
           if(sj != noj) dj++;
          }
          if(si != noi) di++;
     }
     return result;
}
ostream& operator<<(ostream& os, const Matrix& M) {
     for(int i = 0; i < M.rows; i++) {
          for(int j = 0; j < M.cols; j++)
               os << setw(3) << M.m[i][j];
      os << endl;
     }
     return os;
}
// ************************************************
int det(const Matrix& M) {
       int result = 0;
        if(M.rows != M.cols) {
          cout << "Solo matrici quadrate" << endl;
          return 0;
         }
         if(M.cols == 1)
              return M.m[0][0]; // Condizione di stop
         for(int j = 0; j < M.cols; j++) {
               if(j%2)  // colonna dispari
                   result -= M.m[0][j] * det(M.minore[0][j]);
              else if(!(j%2)) // colonna pari
                   result += M.m[0][j] * det(M.minore[0][j]);
         }
            return result;
}
// ************************************************
int main() {
    for(int i = 0; i < 25; i++)
 cout << "\n";
    int dim;
    cout << "Dimensione? ";
                cin >> dim;
    Matrix A(dim, dim);
 A.set_random();
//    cout << A << endl;
 cout << det(A) << endl;
 return 0;
}

Non sono un grande esperto comunque...
generalmente si cercano gli autovalori diagonalizzando e poi si moltiplica.
Per la ricerca degli autovalori ci sono metodi esatti e approssimati.
Generalmente il criterio di convergenza è che la matrice sia strictly
diagonally dominant (scusa ma una traduzione civile non mi viene in mente)
che vuol dire che
abs(a(i,,j)))>sum(a(i,j)))
spero si capisca.
Ad ogni modo visto che qua tiro notte se ti illustro tutti gli algoritmi e
poi faccio su un casino ti do i nomi... con una ricerca su qualche motore
trovi tutto spiegato per bene
Power method (approssimato) miglioramenti per convergenza veloce (Aitken,
Rayleigh)
Deflation (altri autovalori)
Inverse Power method
Jacobi (approssimato??)
Da quello che so io spesso sulle matrici huge si fa affidamento su qualche
caratteristica della matrice per sfruttare le entries vuote, simmetria
etc...
Comunque la cosa non ha UNA maniera per essere risolta quindi ti consiglio
di sbirciare un libro di analisi numerica o di FORTRAN dove generalmente
trovi queste cose
buon divertimento

}
Tenta di lavorare direttamente con i puntatori, piuttosto che con le
matrici: l'algebra dei puntatori è più efficiente.
Infine, tenta di trovare una versione iterativa (se esiste)
dell'algoritmo: laricorsività è molto pesante, nonostante sia più
facile da gestire.
Buon divertimento :-)
_____________________________________________________
Francesco Trotta
http://space.tin.it/computer/ftrotta
elimina "x_" dal mio indirizzo per le risposte

Dimenticavo: per altri algoritmi ne troverai sicuramente ricercando il
libro "Numerical Recipes in C" (non so se esista anche per C++) di cui
esiste una versione in formato pdf scaricabile


abs(a(i,,j)))>sum(a(i,j)))
No...
e`   |a(i,i)| > sum|(a(i,j)|                              per ogni i  
                       i!=j
Comunque se non ricordo male il determinante di una matrice e` dato dal
prodotto dei suoi autovalori, quindi un qualunque metodo (diretto o 
iterativo) per il calcolo degli autovalori puo` essere utilizzato allo
scopo. La regola di Laplace se non erro (non vorrei spararla grossa)
ha complessita` superpolinomiale, mentre ad esempio il metodo
di Gaus e` O((n^3)/3)...
Paolo
-----------------------------------------------------------------
  Paolo Marotta
  Computer Science, University of Pisa
  home:   
http://www.cli.di.unipi.it/~marotta
-----------------------------------------------------------------

 per sapere se questa assunzione sia sempre corretta
Si, dato che vero è qualsiasi valore diverso da zero.
EMG

 di Gaus e` O((n^3)/3)...
Esatto, bravo!
Meno male che ci ha pensato qualcun altro, avevo paura di dover riprendere
in mano io i libri!!!
--
Alessandro Monti

Si, dato che vero č qualsiasi valore diverso da zero.
Siamo d'accordo che ogni valore diverso da 0 sia vero, ma la domanda č
differente. Il valore di un'espressione logica che risulta vera č 1 in
C, sempre? Per intenderci potrebbe esserci il compilatore pluto che
fornisce come valore vero per i confronti il valore 7. Non che sia
molto sensato, ma volendo fare i pignoli...


 Si, dato che vero è qualsiasi valore diverso da zero.
Sbagliato.  Tu stai parlando di tests, tipo "int i; ... if (i)", in cui
l'intero "i" deve essere convertito in bool; ma qui si tratta di una
conversione da bool, (sj != noj), ad int.  Il valore "true" e'
garantito essere convertito in 1, e "false" in 0 (cfr. Stroustrup 3,
pag. 71).
-- 
Maurizio Loreti                      
http://wwwcdf.pd.infn.it/~loreti/mlo.html
Un. of Padova, Dept. of Physics - Padova, Italy          
lor...@padova.infn.it

 molto sensato, ma volendo fare i pignoli...
Ah, ecco. No, non mi sembra si possa fare affidamento sul valore di
un'espressione logica vera. L'importante è che sia diversa da zero e il
compilatore potrebbe produrre il valore che trova più comodo. Forse
sarebbe meglio utilizzare il tipo bool.
EMG

 pag. 71).
Uhm... se č cosě allora questa č un'altra importante differenza tra C e
C++.
EMG

 codice che invece non si anno con le precedenti.
La presenza di salti dipende solo dal set di istruzioni del
microprocessore e dal grado di ottimizzazione del compilatore,
non da qualle delle due istruzioni viene usata.
E non sempre i salti sono piu' lenti, a volte il processore
riesce a sapere abbastanza presto se il salto viene fatto o
no, e a mantenere sempre piena la pipeline.
Un buon compilatore genera esattamente lo stesso codice per
entrambe le istruzioni, un compilatore meno buono ne ottimizza
bene solo una, ma non si puo' sapere quale delle due (a meno
di andare a vedere il codice).
Ho provato a vedere cosa succede con egcs su pentium, compilando
queste due funzioni (ottimizzazione: -O2):
  void f1( int sj, int noj, int & dj )
  {
    if (sj != noj) dj++;
  }
  void f2( int sj, int noj, int & dj )
  {
    dj += (sj != noj);
  }
Il risultato e':
 f1__FiiRi:  // quella con l'if
    movl 8(%esp),%eax         //  eax <- noj
    movl 12(%esp),%edx        //  edx <- &dj
    cmpl %eax,4(%esp)         //  compare( eax, sj )
    je .L434                  //  salta se uguali
    incl (%edx)               //  (*edx)++
 .L434:
    ret
 f2__FiiRi:  // senza l'if, meno leggibile
    movl 12(%esp),%ecx        //  ecx <- &dj
    movl 8(%esp),%eax         //  eax <- noj
    movl (%ecx),%edx          //  edx <- *ecx, cioe' dj
    cmpl %eax,4(%esp)         //  compare( eax, sj )
    je .L436                  //  salto
    incl %edx                 //  edx++ (cioe' dj++)
 .L436:
    movl %edx,(%ecx)          //  dj <- edx
    ret
Il salto e' presente in entrambi i casi, l'unica differenza e'
che nel primo caso dj viene incrementato in memoria, nel secondo
viene prima caricato in un registro e rimesso in memoria alla
fine.  Quando sj==noj (niente incremento), la prima versione
e' piu' veloce, quando sj!=noj i tempi sono gli stessi.
Probabilemte, se invece di funzioni isolate si usano istruzioni
"normali", con molte variabili gia' nei registri, le due
versioni diventano identiche.
Con tutte le ottimizzazioni abilitate, entrambe le funzioni
vengono eliminate e il compilatore mette in dj direttamente
il valore finale.
	Mauro

 C++.
Anche in C un'espressione logica viene valutata come 0 od 1; guarda le
FAQ, Q/A 9.2 .  
Il C++ ha un tipo "bool"; l'operator == (o, se preferisci, il metodo
is_equal) per ogni classe restituisce un bool; e le regole di
conversione da bool ad int sono specificate dal linguaggio.
-- 
Maurizio Loreti                      
http://wwwcdf.pd.infn.it/~loreti/mlo.html
Univ. of Padova, Dept. of Physics - Padova, Italy            
lor...@pd.infn.it


Ragazzi che efficienza. Grazie per i chiarimenti.
Dato che mi sono cadute tutte le ottimizzazioni a livello di codice
che mi erano state riportate, quale strada devo seguire nel creare un
programma?
Tralasciando casi patologici di variabili e codice palesemente
ridondanti, posso supporre che i compilatori abbiano raggiunto un
livello tale da ottimizzare anche porzioni di codice più estese?
Posso quindi tralasciare pseudo-giochetti che non farebbero altro che
rendere il codice meno leggibile (vedi sopra)?
A questo punto verrebbero a cadere strutture di programmazione come
l'espressione condizionale (che può essere tranquillamente convertita
in un costrutto if)?


 FAQ, Q/A 9.2 .
E' vero. Chiedo venia. Mi sono confuso.
EMG

 in un costrutto if)?
L'esistenza dell'operatore ternario non mi pare sia mai stata
giustificata per la maggiore efficienza rispetto all'if ma per la
maggiore concisione. Oggi sappiamo che una delle cose più importanti in
un sorgente è la sua cristallina leggibilità.
EMG

analisi numerica, ma ho un idiosincrasia per il FORTRAN.
Beh... mi pare sia nato (fortran) apposta per quello!


in un costrutto if)?
IMO ha il vantaggio di essere piu` concisa dell'if e data la sua 
unicita` (unico operatore ternario) e l'impossibilita` di eseguirne
l'overloading, puo` essere ancora una scelta conveniente rispetto
all'if (basta non mescolarla in espressioni piu` complesse).


scopo. La regola di Laplace se non erro (non vorrei spararla
Era quello che intendevo dire...
BTW hai un buon libro di analisi numerica da suggerire?
Non mi hanno parlato molto bene di Numerical Recipes.
Pare che ogni tanto ci siano algoritmi non ottimali presentati come tali, e
visto che per il momento mi sarei accontentato di scopiazzare, ed
eventualmente ma solo incidentalmente di capire ;), una roba tipo recipes
andava bene.
I libri di FORTRAN hanno la tendenza ad essere farciti di algoritmi di
analisi numerica, ma ho un idiosincrasia per il FORTRAN.
Grazie per gli eventuali suggerimenti.

Laplace stesso). 
Che guarda caso sono gli autovalori...


Non serve che calcoli gli autovalori : basta che triangolarizzi la matrice e
moltiplichi gli elementi della diagonale (e' un'applicazione del metodo di
Laplace stesso). Per diagonalizzare basta il motodo di Gauss.
Sul mio sito trovi una classe matrice che (tra le altre cose) calcola il det
in questo modo.
Tieni presente che esistono metodi piu' veloci, ma piu' complessi che hanno
una complessitą piu' bassa.
Andrea Murru
http://www.pegacity.it/informatica/case/473/index.htm

...
