Accueil > > > CLASS MATRICE AVEC TEMPLATE
CLASS MATRICE AVEC TEMPLATE
Information sur la source
Description
Cette class matrice réalisé avec les template permet de définir une matrice et d'effectuer la plupart des opérations classique. quelques méthodes sont définies tel que calcul du déterminant, de l'inverse de la trace. J'ai tenté une gestion des exceptions, mais cela est mal fait. J'ai aussi tenté une diagonalisation, mais je me suis heurter à des difficultés d'ordre mathématique (extraction de racines de polynomes) Toute remarque pour compléter ou corriger le fichier est la bien venue. La class a été testé assez largement, et je n'ai jamais constaté d'erreur.
Source
- /* LAMOME Julien 19/10/07 14:12
- ; voir notice ci dessous */
-
- /*
- Utilisation de la struct ci dessous.
- Déclaration :
- - MATRIX<type> nom_matrice(n-ligne,m-colonne) crée une matrice de n lignes sur m colonne d'éléments 'type' de valeurs nul
- - MATRIX<type> nom_matrice(n-ligne/colonne) comme la présédente, mais une matrice carré.
- - MATRIX<type> nom_matrice crée une matrice vide
- Acces données :
- lecture/écriture
- donnée de la matrice : nom_matrice(ligne,colonne)
- données de la matrice : nom_matrice[colonne][ligne] ATTENTION, la convention est inverse
- lecture
- sous matrice : nom_matrice(l1,l2,c1,c2) renvoie la matrice allant de la ligne l1 à l2 et de la colonne c1 à c2
- nombre de ligne : nom_matrice.size_M().n
- nombre de colonnes : nom_matrice.size_M().m
- donnée de la matrice : nom_matrice(i,'r'|'c') renvoi la ligne(r) ou colonne(c) n° i
- Fonctions membres :
- Id() : transforme la matrice en matrice identité, si elle est carré, sinon, ne fait rien
- inverse() : calcul l'inverse d'une matrice carré, quitte le programme sinon
- transpose() : prend la transposer de la matrice.
- affiche() : affiche la matrice, soit avec <stdio.h> soit avec <iostream> include précédent
- det() : calcul le déterminant de la matrice si carrée, renvoi -1 sinon.
- size() : retourne le nombre de lignes d'une matrice carré, et le nombre d'éléments d'une non carré
- size_M() : retourne la taille de la matrice dans la structure de type M_taille.
- Tr() : calcul la trace de la matrice (somme diagonal), renvoi 0 si non carré.
- clear() : rend la matrice vide.
- genre() : retourne le type de matrice : IS_VECTOR, si une des dimmension est de 1; IS_SCALAR, si les 2 dim sont à 1;
- VIDE si les dim sont nul; IS_MATRICE ; et MATRICE_DIM_ERR si les dimensions sont indéterminés
- push_back(type) : rempli une matrice par la fin, la rend carré. !! UTILISATION DECONSEILLER !!
- fct(T fonc(T)) : applique la fonction fonc a chaque membre de la matrice
- Supporte :
- - ostream<<nom_matrice
- - ostream<<nom_matrice.size_M()
- - la plupart des opération algébrique de base : +,-,*,=
- - +,* par "type", et par "vector<type>"
- - les fonctions inverse(matrice), transpose(matrice).
- - la "procedure" TLU, qui à partir de la premiere matrice : A, retourne les matrice triangulaires basse et haute, L, et U, tel que A=L*U
- Compilation tester avec succes avec gcc, jusqu'à la version 3.3.1 compatible cygwin */
-
- #ifndef _MATRIX
- #define _MATRIX
-
- #include <vector>
- #include <cmath>
- #ifdef _STDIO_H_
- #define __ENVOI(x) printf(x);
- #define __ENVOI2(x,y) printf(x,y);
- #else
- #if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
- #define __ENVOI(x) std::cout<<x;
- #else
- #define __ENVOI(x) throw (x);
- #endif
- #endif
- struct M_taille
- {
- int n;
- int m;
- };
- #define IS_MATRICE 2
- #define IS_VECTOR 1
- #define IS_SCALAR 0
- #define VIDE 3
- #define MATRICE_DIM_ERR -1
- template <class T>struct MATRIX
- {
- // variable de la matrice : ses données + sa taille
- private:
- int n,m;
- std::vector<T> M;
- // constructeur, destructeur, operateur, et fonctions
- public:
- MATRIX() {n=0;m=0;std::vector<T> M;};
- MATRIX(int ,int =0); //constructeur
- MATRIX(std::vector<T>); //transformation d'un vecteur en matrice (1,n)
- MATRIX(T);// Attention si T=int !!!!! mettre "MATRIX operator=(T)" ne marche pas
- ~MATRIX(){M.clear();M.reserve(0);};
- int size() const;
- M_taille size_M() const;
- void clear();
- T& operator()(int ,int );
- const T operator()(int ,int )const;
- MATRIX operator()(int ,int ,int ,int );
- std::vector<T> operator()(int,char);
- // std::vector<T> operator[](int);
- T* operator[](int);
- const T* operator[](int)const;
- MATRIX operator=(MATRIX );
- MATRIX operator+=(MATRIX );
- MATRIX operator-=(MATRIX );
- MATRIX operator*=(MATRIX );
- MATRIX operator-() const ;
- MATRIX operator+() const {return *this;};
- MATRIX<T> fct(T t(T)) const;// pour appliquer n'importe quelle fonction
- T Tr() const;
- MATRIX& push_back(T);
- int genre() const;
- T det() const;
- MATRIX& Id();
- MATRIX& inverse();
- MATRIX& transpose();
- void affiche() const;
- void resize(int,int);
- // Fonctions amies
- template <class _T>friend MATRIX<_T> operator+(const MATRIX<_T>&,_T);
- template <class _T>friend MATRIX<_T> operator*(const MATRIX<_T>&,const _T);
- template <class _T>friend MATRIX<_T> operator*(MATRIX<_T>&,MATRIX<_T>&);
- template <class _T>friend MATRIX<_T> operator+(MATRIX<_T>&,MATRIX<_T>&);
- template <class _T>friend MATRIX<_T> operator&(MATRIX<_T>&,MATRIX<_T>&);
- template <class _T>friend MATRIX<_T> operator-(MATRIX<_T>&,MATRIX<_T>&);
- template <class _T>friend bool operator==(MATRIX<_T>,MATRIX<_T>);
- template <class _T>friend std::vector<_T> operator*(MATRIX<_T>,std::vector<_T>);
- template <class _T>friend int TLU(const MATRIX<_T>&,MATRIX<_T>&,MATRIX<_T>&);
- // template <class _T>friend std::ostream& operator<<(std::ostream&, MATRIX<_T>);
- // template <class _T>friend MATRIX<_T>& operator>>(std::istream&, MATRIX<_T>&);
- template <class _T,typename T1>friend MATRIX<_T>& operator>>(T1&, MATRIX<_T>&);
- template <class T1,typename _T>friend T1& operator<<(T1&, MATRIX<_T>);
- };
- /*------------fin de la struct MATRIX-----------------------------------------*/
- /*------------------opération binaire sur les matrices------------------------*/
- template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const T );
- template <class T>MATRIX<T> operator+(const T ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const T );
- template <class T>MATRIX<T> operator*(const T ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator&(MATRIX<T>& ,MATRIX<T>& ); //multiplication terme à terme
- template <class T>MATRIX<T> operator&(const MATRIX<T>& ,const MATRIX<T>& );
- template <class T>MATRIX<T> operator*(MATRIX<T>& ,MATRIX<T>& );
- template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const MATRIX<T>& );
- template<class T> bool operator==(const MATRIX<T>,const MATRIX<T>);
- template<class T> bool operator!=(MATRIX<T>,MATRIX<T>);
- /*------------------opération binaire matrice vecteur------------------------*/
- template <class T>std::vector<T> operator*(MATRIX<T> ,std::vector<T> );
- template <class T>std::vector<T> operator*(std::vector<T> ,MATRIX<T> );
- /*-----------fonction supplémentaire utile au calculs de matrices-------------*/
- template<class T>std::vector<T> operator+(std::vector<T> ,std::vector<T>);
- template<class T>std::vector<T> operator*(T ,std::vector<T>);
- template<class T>std::vector<T> operator*(std::vector<T>,T );
- template<class T>MATRIX<T> transpose(MATRIX<T> );
- template<class T>MATRIX<T> inverse(MATRIX<T> );
- template<class T> int TLU(const MATRIX<T>& A,MATRIX<T>& L,MATRIX<T>& U);
- #if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
- //template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A);
- //std::ostream& operator<<(std::ostream& s, M_taille A);
- #endif // IOSTREAM
- /*--------------------- fonctions de flux pour matrices ----------------------*/
- template <class _T,typename T1> MATRIX<_T>& operator>>(T1&, MATRIX<_T>&);
- template <class T1,typename _T> T1& operator<<(T1&, MATRIX<_T>);
-
-
-
- typedef MATRIX<double> MATRIX2;
- /*----------------------------------------------------------------------------*/
- //--------début des implementation--------------------------------------------*/
- /*----------------------------------------------------------------------------*/
- template<class T>void MATRIX<T>::affiche() const
- {
- #ifdef _STDIO_H_
- for (int i=0;i<n;i++)
- {
- printf("| ");
- for (int j=0;j<m;j++)
- printf("%.1e | ",M[i+n*j]);
- printf("\n");
- }
- #else
- #if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
- for (int i=0;i<n;i++)
- {
- std::cout<<"| ";
- for (int j=0;j<m;j++)
- std::cout<<M[i+n*j]<<" | ";
- std::cout<<"\n";
- }
- #endif
- #endif
- }
- /*------------------opération binaire sur les matrices------------------------*/
- template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const T B)
- {
- MATRIX<T> tmp;
- if (A.m==A.n)
- {
- for (int i=0;i<A.M.size();i++)
- tmp.push_back(A.M[i]+B);
- return tmp;
- }
- MATRIX<T> tmp1(A.n,A.m);
- for (int j=0;j<A.m;j++)
- for (int i=0;i<A.n;i++)
- tmp1(i,j)=A(i,j)+B;
- return tmp1;
- }
- template <class T>MATRIX<T> operator+(const T B,const MATRIX<T>& A)
- {
- return A+B;
- }
- template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const T B)
- {
- MATRIX<T> tmp;
- if (A.m==A.n)
- {
- for (unsigned int i=0;i<A.M.size();i++)
- tmp.push_back(A.M[i]*B);
- return tmp;
- }
- MATRIX<T> tmp1(A.n,A.m);
- for (int j=0;j<A.m;j++)
- for (int i=0;i<A.n;i++)
- tmp1(i,j)=A(i,j)*B;
- return tmp1;
- }
- template <class T>MATRIX<T> operator*(const T B,const MATRIX<T>& A)
- {
- return A*B;
- }
- template <class T>MATRIX<T> operator+(MATRIX<T>& A,MATRIX<T>& B)
- {
- if (A.m==B.m & A.n==B.n) {
- MATRIX<T> tmp;
- for (int i=0;i<A.M.size();i++)
- tmp.push_back(A.M[i]+B.M[i]);
- tmp.n=A.n;
- tmp.m=A.m;
- return tmp; }
- __ENVOI("problème de dimension (operator+)");
- return -1;
- }
- template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const MATRIX<T>& B)
- {
- MATRIX<T> a=A;
- MATRIX<T> b=B;
- return a+b;
- }
- template <class T>MATRIX<T> operator-(MATRIX<T>& A,MATRIX<T>& B)
- {
- if (A.m==B.m & A.n==B.n) {
- MATRIX<T> tmp;
- for (int i=0;i<A.M.size();i++)
- tmp.push_back(A.M[i]-B.M[i]);
- tmp.n=A.n;
- tmp.m=A.m;
- return tmp; }
- __ENVOI("problème de dimension (operator-)");
- return -1;
- }
- template <class T>MATRIX<T> operator-(const MATRIX<T>& A,const MATRIX<T>& B)
- {
- MATRIX<T> a=A;
- MATRIX<T> b=B;
- return a-b;
- }
- template <class T>MATRIX<T> operator&(MATRIX<T>& A,MATRIX<T>& B)
- {
- if (A.m==B.m & A.n==B.n) {
- MATRIX<T> tmp;
- T a;
- for (int i=0;i<A.n*A.m;i++)
- {
- a=A.M[i]*B.M[i];
- tmp.push_back(a);
- }
- tmp.n=A.n;
- tmp.m=A.m;
- return tmp; }
- return -1;
- }
- template <class T>MATRIX<T> operator&(const MATRIX<T>& A,const MATRIX<T>& B)
- {
- MATRIX<T> a=A;
- MATRIX<T> b=B;
- return a&b;
- }
- template <class T>MATRIX<T> operator*(MATRIX<T>& A,MATRIX<T>& B)
- {
- if (A.m==B.n)
- {
- MATRIX<T> tmp;
- T a=0;
- for (int j=0;j<B.m;j++)
- for (int i=0;i<A.n;i++)
- {
- for (int k=0;k<A.m;k++)
- a=a+A(i,k)*B(k,j);
- tmp.push_back(a);
- a=0;
- }
- tmp.n=A.n;
- tmp.m=B.m;
- return (tmp);
- }
- if ((A.n==B.m)&(A.m==1)&(B.n==1)) //très con, reviens au test précedent Am=1 et Bn=1 => Am==Bn !
- {
- MATRIX<T> tmp;
- for (int i=0;i<A.n;i++)
- for (int j=0;j<B.m;j++)
- tmp.push_back(A(i,0)*B(0,j));
- return tmp;
- }
- if (A.n==1 & A.m==1)
- return A(0,0)*B;
- if (B.n==1 & B.m==1)
- return B*A;
- __ENVOI("multiplication de matrices de dimension incompatible\n");
- return -1;
- }
- template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const MATRIX<T>& B)
- {
- MATRIX<T> a=A;
- MATRIX<T> b=B;
- return a*b;
- }
- template<class T> bool operator==(MATRIX<T> A,MATRIX<T> B)
- {
- if( (A.n!=B.n) | (A.m!=B.m) ) return false;
- for (int i=0;i<A.n*A.m;i++)
- if(A.M[i]!=B.M[i]) return false;
- return true;
- }
- template<class T> bool operator!=(MATRIX<T> A,MATRIX<T> B){return !(A==B);}
- /*------------------opération binaire matrice vecteur------------------------*/
- template <class T>std::vector<T> operator*(MATRIX<T> A,std::vector<T> v)
- {
- if (v.size()!=A.n)
- {
- __ENVOI("vecteur de mauvaise dimension\n");
- }
- std::vector<T> tmp;
- T a ;
- for (int i=0;i<A.n;i++)
- {
- a=0;
- for (int j=0;j<A.m;j++)
- a=a+A(i,j)*v[j];
- tmp.push_back(a);
- }
- return tmp;
- }
- template <class T>std::vector<T> operator*(std::vector<T> v,MATRIX<T> A)
- {
- return A*v;
- }
- /*---définitions des fonctions et constructeurs de la strut MATRIX------------*/
- template <class T>MATRIX<T>::MATRIX<T>(int a,int b)
- {
- n=a;
- m=b;
- if (b==0) m=n;
- for (int i=0;i<n*m;i++)
- M.push_back(T(0));
- }
- template <class T>MATRIX<T>::MATRIX<T>(std::vector<T> v)
- {
- n=v.size();
- m=1;
- M=v;
- }
- template <class T>MATRIX<T>::MATRIX<T>(T v)
- {
- n=1;
- m=1;
- M.push_back(v);
- }
- template<class T> T& MATRIX<T>::operator()(int x,int y)
- {
- if ( (x<n)&(y<m) )
- return M[x+n*y];
- __ENVOI("indice trop grand [operator()]\n");
- exit(-1);
- }
- template<class T> const T MATRIX<T>::operator()(int x,int y)const
- {
- if ( (x<n)&(y<m) )
- return M[x+n*y];
- __ENVOI("indice trop grand [operator()]\n");
- exit(-1);
- }
- template<class T> MATRIX<T> MATRIX<T>::operator()(int a,int b,int c,int d)
- {
- if (a>b | b>n | c>d | d>m) {__ENVOI("erreur indice sous matrice");return 0;}
- MATRIX<T> tmp(b-a+1,d-c+1);
- for (int i=0;i<b-a+1;i++)
- for (int j=0;j<d-c+1;j++)
- tmp(i,j)=M[i+a+n*(j+c)];
- return tmp;
- }
- template<class T> std::vector<T> MATRIX<T>::operator()(int x,char c)
- {
- std::vector<T> tmp;
- if (c=='c')
- {
- for (int i=0;i<n;i++)
- tmp.push_back(M[i+n*x]);
- return tmp;
- }
- if (c=='r')
- {
- for (int i=0;i<m;i++)
- tmp.push_back(M[x+n*i]);
- return tmp;
- }
- return tmp;
- }
- /*template<class T> std::vector<T> MATRIX<T>::operator[](int r)
- {
- std::vector<T> tmp(m);
- for (int i=0;i<m;i++)
- tmp[i]=M[i+n*r];
- return tmp;
- }*/
- template<class T> T* MATRIX<T>::operator[](int r) //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
- {
- int nr=n*r;
- return &M[nr];
- }
- template<class T>const T* MATRIX<T>::operator[](int r)const //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
- {
- int nr=n*r;
- return &M[nr];
- }
- template <class T> MATRIX<T> MATRIX<T>::operator=(MATRIX S)
- {
- M=S.M;
- n=S.n;
- m=S.m;
- return *this;
- }
- template <class T> MATRIX<T> MATRIX<T>::operator+=(MATRIX S){*this=*this+S;return *this;}
- template <class T> MATRIX<T> MATRIX<T>::operator-=(MATRIX S){*this=*this-S;return *this;}
- template <class T> MATRIX<T> MATRIX<T>::operator*=(MATRIX S){*this=*this*S;return *this;}
- template <class T> MATRIX<T> MATRIX<T>::operator-()const {return T(-1)**this;}
- template <class T> int MATRIX<T>::size() const
- {
- if (M.size()!=n*m) return -1;
- if (m==n) return n;
- return n*m;
- }
- template <class T>MATRIX<T>& MATRIX<T>::Id()
- {
- if (m==n)
- for (int i=0;i<n;i++)
- M[i+n*i]=1.0;
- return *this;
- }
- template <class T> void MATRIX<T>::clear()
- {
- M.clear();
- n=0;
- m=0;
- }
- template <class T> T MATRIX<T>::Tr() const
- {
- T tr=0;
- for (int i=0;i<n;i++)
- tr+=M[i*(1+n)];
- if (m==n)
- return tr;
- return 0;
- }
- template <class T>MATRIX<T>& MATRIX<T>::inverse()
- {
- if (m!=n)
- {
- __ENVOI("inversion de matrices non carré non implémenter!\n");
- return *this;
- }
- MATRIX<T> MM(n),M1(n),Mi1(n),Minv(n);
- MM=*this;
- M1=MM;
- Mi1.Id();
- Minv=Mi1;
- for (int i=0;i<n;i++)
- {
- for (int j=0;j<n;j++)
- {
- if (MM(i,i)==0)
- {
- __ENVOI("inversion pivot GAUSS impossible : division par 0\n");
- return *this;
- }
- M1(i,j)=MM(i,j)/MM(i,i);
- Mi1(i,j)=Minv(i,j)/(MM(i,i));
- }
- MM=M1;
- Minv=Mi1;
- for (int k=0;k<n;k++)// mise zéros de la colonne
- if (k!=i)
- for (int j=0;j<n;j++)
- {
- M1(k,j)=MM(k,j)-MM(i,j)*MM(k,i);
- Mi1(k,j)=Minv(k,j)-Minv(i,j)*MM(k,i);
- }
- MM=M1;
- Minv=Mi1;
- /*for (int k=0;k<n;k++)//mise à zéro de la ligne. ap essai : pas la peine
- if (k!=i)
- for (int j=0;j<n;j++)
- {
- M1(j,k)=MM(j,k)-MM(i,j)*MM(i,k);
- Mi1(j,k)=Minv(j,k)-Minv(i,j)*MM(i,k);
- }
- MM=M1;
- Minv=Mi1;*/
- }
- *this=Minv;
- return *this;
- }
- template <class T> T MATRIX<T>::det() const
- {
- if (m!=n)
- {
- __ENVOI("déterminant de matrices non carré non implémenter!\n");
- return -1;
- }
- if (n==1) return M[0];
- //if (n==2) return M[0]*M[3]-M[1]*M[2];
- T d=0.0;
- MATRIX tmp;
- for (int i=0;i<n;i++)
- {
- for (int j=n;j<n*n;j++)
- if (int((j-i)*1.0/n)!=(j-i)*1.0/n)
- tmp.push_back(M[j]);
- d=d+pow(-1,double(i))*M[i]*tmp.det();
- tmp.clear();
- }
- return d;
- }
- template <class T> MATRIX<T>& MATRIX<T>::transpose()
- {
- MATRIX<T> P(m,n);
- for (int i=0; i<m; i++)
- for (int j=0; j<n; j++)
- { P(i,j)=M[j+n*i];}
- *this=P;
- return *this;
- }
- template <class T> MATRIX<T>& MATRIX<T>::push_back(T x)
- {
- M.push_back(x);
- n=int(sqrt(double(M.size())));
- m=n;
- return *this;
- }
- template <class T> int MATRIX<T>::genre() const
- {
- if ((m==1)&(n==1))
- return IS_SCALAR;
- if ((m==1)|(n==1))
- return IS_VECTOR;
- if ((m==0)&(n==0))
- return VIDE;
- if ((m>1)&(n>1))
- return IS_MATRICE;
- else
- return MATRICE_DIM_ERR;
- }
- template<class T>void MATRIX<T>::resize(int l,int c)
- {
- if(l*c>m*n)
- {__ENVOI("changement de taille impossible : pas assez d'éléments\n");return;}
- if(l*c==m*n)
- {n=l;m=c;return;}
- for(int i=0;i<l*c-m*n;i++)
- M.pop_back();
- }
- template<class T>M_taille MATRIX<T>::size_M() const
- {
- M_taille tmp;
- tmp.n=n;
- tmp.m=m;
- return tmp;
- }
- template <class T>MATRIX<T> MATRIX<T>::fct(T t(T)) const
- {
- MATRIX<T> A=*this;
- for(int i=0;i<m*n;i++)
- A.M[i]=t(A.M[i]);
- return A;
- }
- /*-----------fonction supplémentaire utile au calculs de matrices-------------*/
- template<class T>std::vector<T> operator+(std::vector<T> x,std::vector<T>y)
- {
- std::vector<T> tmp;
- if ( x.size()==y.size() )
- for (int i=0;i<x.size();i++)
- tmp.push_back(x[i]+y[i]);
- return tmp;
- }
- template<class T>std::vector<T> operator-(std::vector<T> x,std::vector<T>y){return x+T(-1.)*y;}
- template<class T>std::vector<T> operator*(T x,std::vector<T>y)
- {
- std::vector<T> tmp;
- for (int i=0;i<y.size();i++)
- tmp.push_back(x*y[i]);
- return tmp;
- }
- template<class T>std::vector<T> operator*(std::vector<T>y,T x)
- {
- return x*y;
- }
- // fonction transpose et inverse fonctionnant par copie volontairement
- // permettant de ne pas modifier A et de les utiliser avec une const MATRIX
- template<class T>MATRIX<T> transpose(MATRIX<T> A)
- {return A.transpose();}
- template<class T>MATRIX<T> inverse(MATRIX<T> A)
- {return A.inverse();}
- template<class T>int TLU(const MATRIX<T>& A,MATRIX<T>& b,MATRIX<T>& c)
- // fonction retournant une matrice triangulaire haute(c), et basse (b)
- // tel que a=b×c
- {
- if (A.n!=A.m) return -1;
- //MATRIX<T> c(A.n),b(A.n);
- b=MATRIX<T>(A.n);c=b;
- T som=0;
- b(0,0)=A(0,0);
- c(0,0)=1;
- for (int j=1;j<A.n;j++)
- {
- c(0,j)=A(0,j)/c(0,0);
- b(j,0)=A(j,0)/b(0,0);
- }
- for (int k=1;k<A.n;k++)
- for (int kk=k;kk<A.n;kk++)
- {
- som=0;
- for (int l=0;l<kk;l++)
- som+=b(kk,l)*c(l,kk);
- b(kk,kk)=A(kk,kk)-som;
- c(kk,kk)=1;
- som=0;
- for (int l=0;l<k;l++)
- som+=b(k,l)*c(l,kk);
- c(k,kk)=(A(k,kk)-som)/b(kk,kk);
- b(kk,k)=(A(kk,k)-som)/c(k,k);
- }
- return 0;
- }
- #if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
- //template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A)
- // {
- // if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
- // for (int i=0;i<A.n;i++)
- // {
- // s<<"\t|";
- // for (int j=0;j<A.m;j++)
- // s<<A(i,j)<<" | ";
- // s<<"\n";
- // }
- // return s;
- // }
- //template<class T>MATRIX<T>& operator>>(std::istream& s, MATRIX<T>& A)
- // {
- // std::cout<<"Saisi de matrice:\nnombre de lignes\n";
- // s>>A.n;
- // std::cout<<"nombre de colonnes\n";
- // s>>A.m;
- // A=MATRIX<T>(A.n,A.m);
- // std::cout<<"Remplisage par ligne puis par colonne\n";
- // for (int i=0;i<A.n*A.m;i++)
- // s>>A.M[i];
- // //A.transpose();
- // return A;
- // }
- std::ostream& operator<<(std::ostream& s, M_taille A)
- {
- s<<"| "<<A.n<<" | "<<A.m<<" | ";
- return s;
- }
- #endif // IOSTREAM
- template <class T,typename T1> MATRIX<T>& operator>>(T1& s, MATRIX<T>& A)
- {
- __ENVOI("Saisi de matrice:\nnombre de lignes\n");
- s>>A.n;
- __ENVOI("nombre de colonnes\n");
- s>>A.m;
- A=MATRIX<T>(A.n,A.m);
- __ENVOI("Remplisage par ligne puis par colonne\n");
- for (int i=0;i<A.n*A.m;i++)
- s>>A.M[i];
- //A.transpose();
- return A;
- }
- template <class T1,typename _T> T1& operator<<(T1& s, MATRIX<_T> A)
- {
- if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
- for (int i=0;i<A.n;i++)
- {
- s<<"\t|";
- for (int j=0;j<A.m;j++)
- s<<A(i,j)<<" | ";
- s<<"\n";
- }
- return s;
- }
- template<class T>MATRIX<T> pow(const MATRIX<T>& r,const MATRIX<T>& e)// à généraliser pour les autres fonctions : sin cos exp ln etc...
- {
- MATRIX<T> tmp1=r;
- if ((r.genre()==IS_SCALAR) & (e.genre()==IS_SCALAR))
- return pow(r(0,0),e(0,0));
- else
- if(e.genre()==IS_SCALAR)
- if(e(0,0)>0)
- for(int i=1;i<e(0,0);i++)
- tmp1=tmp1*r;
- else
- tmp1=inverse(pow(r,T(-1.0)*e));
- return tmp1;
- }
- template<class T>MATRIX<T> cos(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
- {
- MATRIX<T> tmp1=r;
- int n=tmp1.size_M().n,m=tmp1.size_M().m;
- for(int i=0;i<n;i++)
- for(int j=0;j<m;j++)
- tmp1(i,j)=cos(tmp1(i,j));
- return tmp1;
- }
- template<class T>MATRIX<T> log(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
- {
- MATRIX<T> tmp1=r;
- int n=tmp1.size_M().n,m=tmp1.size_M().m;
- for(int i=0;i<n;i++)
- for(int j=0;j<m;j++)
- tmp1(i,j)=log(tmp1(i,j));
- return tmp1;
- }
- template<class T>MATRIX<T> fabs(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
- {
- MATRIX<T> tmp1=r;
- int n=tmp1.size_M().n,m=tmp1.size_M().m;
- for(int i=0;i<n;i++)
- for(int j=0;j<m;j++)
- tmp1(i,j)=fabs(tmp1(i,j));
- return tmp1;
- }
- template<class T>MATRIX<T> acos(const MATRIX<T>& a){return a.fct(acos);}
- template<class T>MATRIX<T> operator<(MATRIX<T> r,MATRIX<T> s)
- // Compare les matrices. Tout d'abord le nombre d'éléments, si s a plus délément s que r
- // Ensuite compare leur déterminant si elles sont carrées
- // Pour finir leur ordre, ainsi 2×2>1×4
- // le cas 1×4<4×1 renvoi zero car on ne peut comparer
- {
- MATRIX<T> tmp1(1);
- if(r.size()<s.size())
- {tmp1(0,0)=1;return tmp1;}
- else if(r.size()>s.size())
- {tmp1(0,0)=0;return tmp1;}
- else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
- {tmp1(0,0)=(T)(r.det()<s.det());return tmp1;}
- int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
- if(fmin(rn,rm)<fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
- tmp1(0,0)=0;
- return tmp1;
- }
- template<class T>MATRIX<T> operator>(MATRIX<T> r,MATRIX<T> s)
- {
- MATRIX<T> tmp1(1);
- if(r.size()>s.size())
- {tmp1(0,0)=1;return tmp1;}
- else if(r.size()<s.size())
- {tmp1(0,0)=0;return tmp1;}
- else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
- {tmp1(0,0)=(T)(r.det()>s.det());return tmp1;}
- int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
- if(fmin(rn,rm)>fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
- tmp1(0,0)=0;
- return tmp1;
- }
- //template<class T>MATRIX<T> operator==(MATRIX<T> r,MATRIX<T> s) //définie par défaut
- // {
- // MATRIX<T> tmp1(1);
- // int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
- // if((rn!=sn)|(rm!=sm)){tmp1(0,0)=0;return tmp1;}
- // else
- // for(int i=0;i<rn;i++)
- // for(int j=0;j<rm;j++)
- // if(r(i,j)!=s(i,j)){tmp1(0,0)=0;return tmp1;}
- // tmp(0,0)=1;
- // return tmp1;
- // }
- #endif //ndef _matrice
- /*int diagonalisation(MATRIX2,MATRIX2&,MATRIX2&,std::vector<double>&);
- int diagonalisation(MATRIX2 A,MATRIX2& D,MATRIX2& P,std::vector<double>& Lambda)
- {
- double s,df;
- D=MATRIX2(A.n);P=MATRIX2(A.n);Lambda.clear();
- if (A.n!=A.m) {__ENVOI("imposible de diagonalise une matrice non carre!\n");exit(0);}
- double inter=0,tmp1,tmp2,tmp3;
- double pr=200;//precision
- for (int i=0;i<A.n*A.n;i++)inter+=fabs(A.M[i]);
- for (double l=-inter;l<inter;l+=inter/pr)
- / *{
- while(1)
- {
- df=((A-l*MATRIX2(A.n).Id()).det()-(A-(l-inter/(20*pr))*MATRIX2(A.n).Id()).det())/(inter/(20*pr));
- s=l+(A-l*MATRIX2(A.n).Id()).det()/df;cout<<s-l<<endl;if(s<l) {l+=inter;continue;}
- if (s-l!=0) l=s;
- else break;if (!-inter<l) system("pause");
- }
- Lambda.push_back(l);
- }//méthode de Newton(proto necessite d'avoir la dérivée du polynome caractéristique)* /
- / * {
- tmp1=(A-l*MATRIX2(A.n).Id()).det();
- if (tmp1*tmp2<0 & l!=-inter)
- Lambda.push_back(l-inter*2.0/pr);
- if (tmp1==0)
- Lambda.push_back(l);
- tmp2=tmp1;
- } //méthode perso, lente, possibilité de sauter des solutions, et impresice
- if (Lambda.size()!=A.n) {cout<<MATRIX2(Lambda);__ENVOI("probleme de nombre de VP\n");exit(0);}
- for (int i=0;i<A.n;i++)D(i,i)=Lambda[i];
- tmp2=0;for (int i=1;i<A.n;i++) {if (A(0,i)!=0)tmp2=1;}tmp3=tmp2;//traitement du cas où la 1ére coordonnée du VP est nulle.(proto)
- MATRIX2 AV(A.n-1),AVt;
- for(int i=0;i<A.n;i++)
- {
- for (int k=0;k<A.n-1;k++)for (int j=0;j<A.n-1;j++)AV(k,j)=A(k+1,j+1);
- AV=AV-Lambda[i]*MATRIX2(AV.n).Id();
- if(abs(A(0,0)-Lambda[i])<2*inter/pr) tmp3=1;else tmp3=tmp2;
- tmp1=AV.det();
- P(0,i)=tmp3;
- for(int j=1;j<A.n;j++)
- {
- AVt=AV;
- for (int k=0;k<AV.n;k++) AVt(k,j-1)=-A(k+1,0)/ **tmp3* /;
- / * P(j,i)=AVt.det()/tmp1;
- }
- }
- return 0;
- }* /
- / * pour passer n m et M en private, il faut definir les fonction qui les utilise en friend dans la declaration
- (si j'ai bien compris). et rajouter eventuellement quelques fonctions, comme resize pour aller avec la fonction
- push_back, si l'on ne veut pas une matrice carré */
/* LAMOME Julien 19/10/07 14:12
; voir notice ci dessous */
/*
Utilisation de la struct ci dessous.
Déclaration :
- MATRIX<type> nom_matrice(n-ligne,m-colonne) crée une matrice de n lignes sur m colonne d'éléments 'type' de valeurs nul
- MATRIX<type> nom_matrice(n-ligne/colonne) comme la présédente, mais une matrice carré.
- MATRIX<type> nom_matrice crée une matrice vide
Acces données :
lecture/écriture
donnée de la matrice : nom_matrice(ligne,colonne)
données de la matrice : nom_matrice[colonne][ligne] ATTENTION, la convention est inverse
lecture
sous matrice : nom_matrice(l1,l2,c1,c2) renvoie la matrice allant de la ligne l1 à l2 et de la colonne c1 à c2
nombre de ligne : nom_matrice.size_M().n
nombre de colonnes : nom_matrice.size_M().m
donnée de la matrice : nom_matrice(i,'r'|'c') renvoi la ligne(r) ou colonne(c) n° i
Fonctions membres :
Id() : transforme la matrice en matrice identité, si elle est carré, sinon, ne fait rien
inverse() : calcul l'inverse d'une matrice carré, quitte le programme sinon
transpose() : prend la transposer de la matrice.
affiche() : affiche la matrice, soit avec <stdio.h> soit avec <iostream> include précédent
det() : calcul le déterminant de la matrice si carrée, renvoi -1 sinon.
size() : retourne le nombre de lignes d'une matrice carré, et le nombre d'éléments d'une non carré
size_M() : retourne la taille de la matrice dans la structure de type M_taille.
Tr() : calcul la trace de la matrice (somme diagonal), renvoi 0 si non carré.
clear() : rend la matrice vide.
genre() : retourne le type de matrice : IS_VECTOR, si une des dimmension est de 1; IS_SCALAR, si les 2 dim sont à 1;
VIDE si les dim sont nul; IS_MATRICE ; et MATRICE_DIM_ERR si les dimensions sont indéterminés
push_back(type) : rempli une matrice par la fin, la rend carré. !! UTILISATION DECONSEILLER !!
fct(T fonc(T)) : applique la fonction fonc a chaque membre de la matrice
Supporte :
- ostream<<nom_matrice
- ostream<<nom_matrice.size_M()
- la plupart des opération algébrique de base : +,-,*,=
- +,* par "type", et par "vector<type>"
- les fonctions inverse(matrice), transpose(matrice).
- la "procedure" TLU, qui à partir de la premiere matrice : A, retourne les matrice triangulaires basse et haute, L, et U, tel que A=L*U
Compilation tester avec succes avec gcc, jusqu'à la version 3.3.1 compatible cygwin */
#ifndef _MATRIX
#define _MATRIX
#include <vector>
#include <cmath>
#ifdef _STDIO_H_
#define __ENVOI(x) printf(x);
#define __ENVOI2(x,y) printf(x,y);
#else
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
#define __ENVOI(x) std::cout<<x;
#else
#define __ENVOI(x) throw (x);
#endif
#endif
struct M_taille
{
int n;
int m;
};
#define IS_MATRICE 2
#define IS_VECTOR 1
#define IS_SCALAR 0
#define VIDE 3
#define MATRICE_DIM_ERR -1
template <class T>struct MATRIX
{
// variable de la matrice : ses données + sa taille
private:
int n,m;
std::vector<T> M;
// constructeur, destructeur, operateur, et fonctions
public:
MATRIX() {n=0;m=0;std::vector<T> M;};
MATRIX(int ,int =0); //constructeur
MATRIX(std::vector<T>); //transformation d'un vecteur en matrice (1,n)
MATRIX(T);// Attention si T=int !!!!! mettre "MATRIX operator=(T)" ne marche pas
~MATRIX(){M.clear();M.reserve(0);};
int size() const;
M_taille size_M() const;
void clear();
T& operator()(int ,int );
const T operator()(int ,int )const;
MATRIX operator()(int ,int ,int ,int );
std::vector<T> operator()(int,char);
// std::vector<T> operator[](int);
T* operator[](int);
const T* operator[](int)const;
MATRIX operator=(MATRIX );
MATRIX operator+=(MATRIX );
MATRIX operator-=(MATRIX );
MATRIX operator*=(MATRIX );
MATRIX operator-() const ;
MATRIX operator+() const {return *this;};
MATRIX<T> fct(T t(T)) const;// pour appliquer n'importe quelle fonction
T Tr() const;
MATRIX& push_back(T);
int genre() const;
T det() const;
MATRIX& Id();
MATRIX& inverse();
MATRIX& transpose();
void affiche() const;
void resize(int,int);
// Fonctions amies
template <class _T>friend MATRIX<_T> operator+(const MATRIX<_T>&,_T);
template <class _T>friend MATRIX<_T> operator*(const MATRIX<_T>&,const _T);
template <class _T>friend MATRIX<_T> operator*(MATRIX<_T>&,MATRIX<_T>&);
template <class _T>friend MATRIX<_T> operator+(MATRIX<_T>&,MATRIX<_T>&);
template <class _T>friend MATRIX<_T> operator&(MATRIX<_T>&,MATRIX<_T>&);
template <class _T>friend MATRIX<_T> operator-(MATRIX<_T>&,MATRIX<_T>&);
template <class _T>friend bool operator==(MATRIX<_T>,MATRIX<_T>);
template <class _T>friend std::vector<_T> operator*(MATRIX<_T>,std::vector<_T>);
template <class _T>friend int TLU(const MATRIX<_T>&,MATRIX<_T>&,MATRIX<_T>&);
// template <class _T>friend std::ostream& operator<<(std::ostream&, MATRIX<_T>);
// template <class _T>friend MATRIX<_T>& operator>>(std::istream&, MATRIX<_T>&);
template <class _T,typename T1>friend MATRIX<_T>& operator>>(T1&, MATRIX<_T>&);
template <class T1,typename _T>friend T1& operator<<(T1&, MATRIX<_T>);
};
/*------------fin de la struct MATRIX-----------------------------------------*/
/*------------------opération binaire sur les matrices------------------------*/
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const T );
template <class T>MATRIX<T> operator+(const T ,const MATRIX<T>& );
template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const T );
template <class T>MATRIX<T> operator*(const T ,const MATRIX<T>& );
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator+(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator-(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator&(MATRIX<T>& ,MATRIX<T>& ); //multiplication terme à terme
template <class T>MATRIX<T> operator&(const MATRIX<T>& ,const MATRIX<T>& );
template <class T>MATRIX<T> operator*(MATRIX<T>& ,MATRIX<T>& );
template <class T>MATRIX<T> operator*(const MATRIX<T>& ,const MATRIX<T>& );
template<class T> bool operator==(const MATRIX<T>,const MATRIX<T>);
template<class T> bool operator!=(MATRIX<T>,MATRIX<T>);
/*------------------opération binaire matrice vecteur------------------------*/
template <class T>std::vector<T> operator*(MATRIX<T> ,std::vector<T> );
template <class T>std::vector<T> operator*(std::vector<T> ,MATRIX<T> );
/*-----------fonction supplémentaire utile au calculs de matrices-------------*/
template<class T>std::vector<T> operator+(std::vector<T> ,std::vector<T>);
template<class T>std::vector<T> operator*(T ,std::vector<T>);
template<class T>std::vector<T> operator*(std::vector<T>,T );
template<class T>MATRIX<T> transpose(MATRIX<T> );
template<class T>MATRIX<T> inverse(MATRIX<T> );
template<class T> int TLU(const MATRIX<T>& A,MATRIX<T>& L,MATRIX<T>& U);
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
//template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A);
//std::ostream& operator<<(std::ostream& s, M_taille A);
#endif // IOSTREAM
/*--------------------- fonctions de flux pour matrices ----------------------*/
template <class _T,typename T1> MATRIX<_T>& operator>>(T1&, MATRIX<_T>&);
template <class T1,typename _T> T1& operator<<(T1&, MATRIX<_T>);
typedef MATRIX<double> MATRIX2;
/*----------------------------------------------------------------------------*/
//--------début des implementation--------------------------------------------*/
/*----------------------------------------------------------------------------*/
template<class T>void MATRIX<T>::affiche() const
{
#ifdef _STDIO_H_
for (int i=0;i<n;i++)
{
printf("| ");
for (int j=0;j<m;j++)
printf("%.1e | ",M[i+n*j]);
printf("\n");
}
#else
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
for (int i=0;i<n;i++)
{
std::cout<<"| ";
for (int j=0;j<m;j++)
std::cout<<M[i+n*j]<<" | ";
std::cout<<"\n";
}
#endif
#endif
}
/*------------------opération binaire sur les matrices------------------------*/
template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const T B)
{
MATRIX<T> tmp;
if (A.m==A.n)
{
for (int i=0;i<A.M.size();i++)
tmp.push_back(A.M[i]+B);
return tmp;
}
MATRIX<T> tmp1(A.n,A.m);
for (int j=0;j<A.m;j++)
for (int i=0;i<A.n;i++)
tmp1(i,j)=A(i,j)+B;
return tmp1;
}
template <class T>MATRIX<T> operator+(const T B,const MATRIX<T>& A)
{
return A+B;
}
template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const T B)
{
MATRIX<T> tmp;
if (A.m==A.n)
{
for (unsigned int i=0;i<A.M.size();i++)
tmp.push_back(A.M[i]*B);
return tmp;
}
MATRIX<T> tmp1(A.n,A.m);
for (int j=0;j<A.m;j++)
for (int i=0;i<A.n;i++)
tmp1(i,j)=A(i,j)*B;
return tmp1;
}
template <class T>MATRIX<T> operator*(const T B,const MATRIX<T>& A)
{
return A*B;
}
template <class T>MATRIX<T> operator+(MATRIX<T>& A,MATRIX<T>& B)
{
if (A.m==B.m & A.n==B.n) {
MATRIX<T> tmp;
for (int i=0;i<A.M.size();i++)
tmp.push_back(A.M[i]+B.M[i]);
tmp.n=A.n;
tmp.m=A.m;
return tmp; }
__ENVOI("problème de dimension (operator+)");
return -1;
}
template <class T>MATRIX<T> operator+(const MATRIX<T>& A,const MATRIX<T>& B)
{
MATRIX<T> a=A;
MATRIX<T> b=B;
return a+b;
}
template <class T>MATRIX<T> operator-(MATRIX<T>& A,MATRIX<T>& B)
{
if (A.m==B.m & A.n==B.n) {
MATRIX<T> tmp;
for (int i=0;i<A.M.size();i++)
tmp.push_back(A.M[i]-B.M[i]);
tmp.n=A.n;
tmp.m=A.m;
return tmp; }
__ENVOI("problème de dimension (operator-)");
return -1;
}
template <class T>MATRIX<T> operator-(const MATRIX<T>& A,const MATRIX<T>& B)
{
MATRIX<T> a=A;
MATRIX<T> b=B;
return a-b;
}
template <class T>MATRIX<T> operator&(MATRIX<T>& A,MATRIX<T>& B)
{
if (A.m==B.m & A.n==B.n) {
MATRIX<T> tmp;
T a;
for (int i=0;i<A.n*A.m;i++)
{
a=A.M[i]*B.M[i];
tmp.push_back(a);
}
tmp.n=A.n;
tmp.m=A.m;
return tmp; }
return -1;
}
template <class T>MATRIX<T> operator&(const MATRIX<T>& A,const MATRIX<T>& B)
{
MATRIX<T> a=A;
MATRIX<T> b=B;
return a&b;
}
template <class T>MATRIX<T> operator*(MATRIX<T>& A,MATRIX<T>& B)
{
if (A.m==B.n)
{
MATRIX<T> tmp;
T a=0;
for (int j=0;j<B.m;j++)
for (int i=0;i<A.n;i++)
{
for (int k=0;k<A.m;k++)
a=a+A(i,k)*B(k,j);
tmp.push_back(a);
a=0;
}
tmp.n=A.n;
tmp.m=B.m;
return (tmp);
}
if ((A.n==B.m)&(A.m==1)&(B.n==1)) //très con, reviens au test précedent Am=1 et Bn=1 => Am==Bn !
{
MATRIX<T> tmp;
for (int i=0;i<A.n;i++)
for (int j=0;j<B.m;j++)
tmp.push_back(A(i,0)*B(0,j));
return tmp;
}
if (A.n==1 & A.m==1)
return A(0,0)*B;
if (B.n==1 & B.m==1)
return B*A;
__ENVOI("multiplication de matrices de dimension incompatible\n");
return -1;
}
template <class T>MATRIX<T> operator*(const MATRIX<T>& A,const MATRIX<T>& B)
{
MATRIX<T> a=A;
MATRIX<T> b=B;
return a*b;
}
template<class T> bool operator==(MATRIX<T> A,MATRIX<T> B)
{
if( (A.n!=B.n) | (A.m!=B.m) ) return false;
for (int i=0;i<A.n*A.m;i++)
if(A.M[i]!=B.M[i]) return false;
return true;
}
template<class T> bool operator!=(MATRIX<T> A,MATRIX<T> B){return !(A==B);}
/*------------------opération binaire matrice vecteur------------------------*/
template <class T>std::vector<T> operator*(MATRIX<T> A,std::vector<T> v)
{
if (v.size()!=A.n)
{
__ENVOI("vecteur de mauvaise dimension\n");
}
std::vector<T> tmp;
T a ;
for (int i=0;i<A.n;i++)
{
a=0;
for (int j=0;j<A.m;j++)
a=a+A(i,j)*v[j];
tmp.push_back(a);
}
return tmp;
}
template <class T>std::vector<T> operator*(std::vector<T> v,MATRIX<T> A)
{
return A*v;
}
/*---définitions des fonctions et constructeurs de la strut MATRIX------------*/
template <class T>MATRIX<T>::MATRIX<T>(int a,int b)
{
n=a;
m=b;
if (b==0) m=n;
for (int i=0;i<n*m;i++)
M.push_back(T(0));
}
template <class T>MATRIX<T>::MATRIX<T>(std::vector<T> v)
{
n=v.size();
m=1;
M=v;
}
template <class T>MATRIX<T>::MATRIX<T>(T v)
{
n=1;
m=1;
M.push_back(v);
}
template<class T> T& MATRIX<T>::operator()(int x,int y)
{
if ( (x<n)&(y<m) )
return M[x+n*y];
__ENVOI("indice trop grand [operator()]\n");
exit(-1);
}
template<class T> const T MATRIX<T>::operator()(int x,int y)const
{
if ( (x<n)&(y<m) )
return M[x+n*y];
__ENVOI("indice trop grand [operator()]\n");
exit(-1);
}
template<class T> MATRIX<T> MATRIX<T>::operator()(int a,int b,int c,int d)
{
if (a>b | b>n | c>d | d>m) {__ENVOI("erreur indice sous matrice");return 0;}
MATRIX<T> tmp(b-a+1,d-c+1);
for (int i=0;i<b-a+1;i++)
for (int j=0;j<d-c+1;j++)
tmp(i,j)=M[i+a+n*(j+c)];
return tmp;
}
template<class T> std::vector<T> MATRIX<T>::operator()(int x,char c)
{
std::vector<T> tmp;
if (c=='c')
{
for (int i=0;i<n;i++)
tmp.push_back(M[i+n*x]);
return tmp;
}
if (c=='r')
{
for (int i=0;i<m;i++)
tmp.push_back(M[x+n*i]);
return tmp;
}
return tmp;
}
/*template<class T> std::vector<T> MATRIX<T>::operator[](int r)
{
std::vector<T> tmp(m);
for (int i=0;i<m;i++)
tmp[i]=M[i+n*r];
return tmp;
}*/
template<class T> T* MATRIX<T>::operator[](int r) //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
{
int nr=n*r;
return &M[nr];
}
template<class T>const T* MATRIX<T>::operator[](int r)const //renvoie l'adresse T* du premier element de la colonne r, donc T[l] sera la ligne l de la colonne r
{
int nr=n*r;
return &M[nr];
}
template <class T> MATRIX<T> MATRIX<T>::operator=(MATRIX S)
{
M=S.M;
n=S.n;
m=S.m;
return *this;
}
template <class T> MATRIX<T> MATRIX<T>::operator+=(MATRIX S){*this=*this+S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator-=(MATRIX S){*this=*this-S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator*=(MATRIX S){*this=*this*S;return *this;}
template <class T> MATRIX<T> MATRIX<T>::operator-()const {return T(-1)**this;}
template <class T> int MATRIX<T>::size() const
{
if (M.size()!=n*m) return -1;
if (m==n) return n;
return n*m;
}
template <class T>MATRIX<T>& MATRIX<T>::Id()
{
if (m==n)
for (int i=0;i<n;i++)
M[i+n*i]=1.0;
return *this;
}
template <class T> void MATRIX<T>::clear()
{
M.clear();
n=0;
m=0;
}
template <class T> T MATRIX<T>::Tr() const
{
T tr=0;
for (int i=0;i<n;i++)
tr+=M[i*(1+n)];
if (m==n)
return tr;
return 0;
}
template <class T>MATRIX<T>& MATRIX<T>::inverse()
{
if (m!=n)
{
__ENVOI("inversion de matrices non carré non implémenter!\n");
return *this;
}
MATRIX<T> MM(n),M1(n),Mi1(n),Minv(n);
MM=*this;
M1=MM;
Mi1.Id();
Minv=Mi1;
for (int i=0;i<n;i++)
{
for (int j=0;j<n;j++)
{
if (MM(i,i)==0)
{
__ENVOI("inversion pivot GAUSS impossible : division par 0\n");
return *this;
}
M1(i,j)=MM(i,j)/MM(i,i);
Mi1(i,j)=Minv(i,j)/(MM(i,i));
}
MM=M1;
Minv=Mi1;
for (int k=0;k<n;k++)// mise zéros de la colonne
if (k!=i)
for (int j=0;j<n;j++)
{
M1(k,j)=MM(k,j)-MM(i,j)*MM(k,i);
Mi1(k,j)=Minv(k,j)-Minv(i,j)*MM(k,i);
}
MM=M1;
Minv=Mi1;
/*for (int k=0;k<n;k++)//mise à zéro de la ligne. ap essai : pas la peine
if (k!=i)
for (int j=0;j<n;j++)
{
M1(j,k)=MM(j,k)-MM(i,j)*MM(i,k);
Mi1(j,k)=Minv(j,k)-Minv(i,j)*MM(i,k);
}
MM=M1;
Minv=Mi1;*/
}
*this=Minv;
return *this;
}
template <class T> T MATRIX<T>::det() const
{
if (m!=n)
{
__ENVOI("déterminant de matrices non carré non implémenter!\n");
return -1;
}
if (n==1) return M[0];
//if (n==2) return M[0]*M[3]-M[1]*M[2];
T d=0.0;
MATRIX tmp;
for (int i=0;i<n;i++)
{
for (int j=n;j<n*n;j++)
if (int((j-i)*1.0/n)!=(j-i)*1.0/n)
tmp.push_back(M[j]);
d=d+pow(-1,double(i))*M[i]*tmp.det();
tmp.clear();
}
return d;
}
template <class T> MATRIX<T>& MATRIX<T>::transpose()
{
MATRIX<T> P(m,n);
for (int i=0; i<m; i++)
for (int j=0; j<n; j++)
{ P(i,j)=M[j+n*i];}
*this=P;
return *this;
}
template <class T> MATRIX<T>& MATRIX<T>::push_back(T x)
{
M.push_back(x);
n=int(sqrt(double(M.size())));
m=n;
return *this;
}
template <class T> int MATRIX<T>::genre() const
{
if ((m==1)&(n==1))
return IS_SCALAR;
if ((m==1)|(n==1))
return IS_VECTOR;
if ((m==0)&(n==0))
return VIDE;
if ((m>1)&(n>1))
return IS_MATRICE;
else
return MATRICE_DIM_ERR;
}
template<class T>void MATRIX<T>::resize(int l,int c)
{
if(l*c>m*n)
{__ENVOI("changement de taille impossible : pas assez d'éléments\n");return;}
if(l*c==m*n)
{n=l;m=c;return;}
for(int i=0;i<l*c-m*n;i++)
M.pop_back();
}
template<class T>M_taille MATRIX<T>::size_M() const
{
M_taille tmp;
tmp.n=n;
tmp.m=m;
return tmp;
}
template <class T>MATRIX<T> MATRIX<T>::fct(T t(T)) const
{
MATRIX<T> A=*this;
for(int i=0;i<m*n;i++)
A.M[i]=t(A.M[i]);
return A;
}
/*-----------fonction supplémentaire utile au calculs de matrices-------------*/
template<class T>std::vector<T> operator+(std::vector<T> x,std::vector<T>y)
{
std::vector<T> tmp;
if ( x.size()==y.size() )
for (int i=0;i<x.size();i++)
tmp.push_back(x[i]+y[i]);
return tmp;
}
template<class T>std::vector<T> operator-(std::vector<T> x,std::vector<T>y){return x+T(-1.)*y;}
template<class T>std::vector<T> operator*(T x,std::vector<T>y)
{
std::vector<T> tmp;
for (int i=0;i<y.size();i++)
tmp.push_back(x*y[i]);
return tmp;
}
template<class T>std::vector<T> operator*(std::vector<T>y,T x)
{
return x*y;
}
// fonction transpose et inverse fonctionnant par copie volontairement
// permettant de ne pas modifier A et de les utiliser avec une const MATRIX
template<class T>MATRIX<T> transpose(MATRIX<T> A)
{return A.transpose();}
template<class T>MATRIX<T> inverse(MATRIX<T> A)
{return A.inverse();}
template<class T>int TLU(const MATRIX<T>& A,MATRIX<T>& b,MATRIX<T>& c)
// fonction retournant une matrice triangulaire haute(c), et basse (b)
// tel que a=b×c
{
if (A.n!=A.m) return -1;
//MATRIX<T> c(A.n),b(A.n);
b=MATRIX<T>(A.n);c=b;
T som=0;
b(0,0)=A(0,0);
c(0,0)=1;
for (int j=1;j<A.n;j++)
{
c(0,j)=A(0,j)/c(0,0);
b(j,0)=A(j,0)/b(0,0);
}
for (int k=1;k<A.n;k++)
for (int kk=k;kk<A.n;kk++)
{
som=0;
for (int l=0;l<kk;l++)
som+=b(kk,l)*c(l,kk);
b(kk,kk)=A(kk,kk)-som;
c(kk,kk)=1;
som=0;
for (int l=0;l<k;l++)
som+=b(k,l)*c(l,kk);
c(k,kk)=(A(k,kk)-som)/b(kk,kk);
b(kk,k)=(A(kk,k)-som)/c(k,k);
}
return 0;
}
#if defined(_CPP_IOSTREAM) || defined(__IOSTREAM__)
//template<class T>std::ostream& operator<<(std::ostream& s, MATRIX<T> A)
// {
// if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
// for (int i=0;i<A.n;i++)
// {
// s<<"\t|";
// for (int j=0;j<A.m;j++)
// s<<A(i,j)<<" | ";
// s<<"\n";
// }
// return s;
// }
//template<class T>MATRIX<T>& operator>>(std::istream& s, MATRIX<T>& A)
// {
// std::cout<<"Saisi de matrice:\nnombre de lignes\n";
// s>>A.n;
// std::cout<<"nombre de colonnes\n";
// s>>A.m;
// A=MATRIX<T>(A.n,A.m);
// std::cout<<"Remplisage par ligne puis par colonne\n";
// for (int i=0;i<A.n*A.m;i++)
// s>>A.M[i];
// //A.transpose();
// return A;
// }
std::ostream& operator<<(std::ostream& s, M_taille A)
{
s<<"| "<<A.n<<" | "<<A.m<<" | ";
return s;
}
#endif // IOSTREAM
template <class T,typename T1> MATRIX<T>& operator>>(T1& s, MATRIX<T>& A)
{
__ENVOI("Saisi de matrice:\nnombre de lignes\n");
s>>A.n;
__ENVOI("nombre de colonnes\n");
s>>A.m;
A=MATRIX<T>(A.n,A.m);
__ENVOI("Remplisage par ligne puis par colonne\n");
for (int i=0;i<A.n*A.m;i++)
s>>A.M[i];
//A.transpose();
return A;
}
template <class T1,typename _T> T1& operator<<(T1& s, MATRIX<_T> A)
{
if(A.genre()==IS_SCALAR) {s<<A(0,0);return s;}
for (int i=0;i<A.n;i++)
{
s<<"\t|";
for (int j=0;j<A.m;j++)
s<<A(i,j)<<" | ";
s<<"\n";
}
return s;
}
template<class T>MATRIX<T> pow(const MATRIX<T>& r,const MATRIX<T>& e)// à généraliser pour les autres fonctions : sin cos exp ln etc...
{
MATRIX<T> tmp1=r;
if ((r.genre()==IS_SCALAR) & (e.genre()==IS_SCALAR))
return pow(r(0,0),e(0,0));
else
if(e.genre()==IS_SCALAR)
if(e(0,0)>0)
for(int i=1;i<e(0,0);i++)
tmp1=tmp1*r;
else
tmp1=inverse(pow(r,T(-1.0)*e));
return tmp1;
}
template<class T>MATRIX<T> cos(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
{
MATRIX<T> tmp1=r;
int n=tmp1.size_M().n,m=tmp1.size_M().m;
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
tmp1(i,j)=cos(tmp1(i,j));
return tmp1;
}
template<class T>MATRIX<T> log(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
{
MATRIX<T> tmp1=r;
int n=tmp1.size_M().n,m=tmp1.size_M().m;
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
tmp1(i,j)=log(tmp1(i,j));
return tmp1;
}
template<class T>MATRIX<T> fabs(const MATRIX<T>& r)// à généraliser pour les autres fonctions : sin cos exp ln etc...
{
MATRIX<T> tmp1=r;
int n=tmp1.size_M().n,m=tmp1.size_M().m;
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
tmp1(i,j)=fabs(tmp1(i,j));
return tmp1;
}
template<class T>MATRIX<T> acos(const MATRIX<T>& a){return a.fct(acos);}
template<class T>MATRIX<T> operator<(MATRIX<T> r,MATRIX<T> s)
// Compare les matrices. Tout d'abord le nombre d'éléments, si s a plus délément s que r
// Ensuite compare leur déterminant si elles sont carrées
// Pour finir leur ordre, ainsi 2×2>1×4
// le cas 1×4<4×1 renvoi zero car on ne peut comparer
{
MATRIX<T> tmp1(1);
if(r.size()<s.size())
{tmp1(0,0)=1;return tmp1;}
else if(r.size()>s.size())
{tmp1(0,0)=0;return tmp1;}
else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
{tmp1(0,0)=(T)(r.det()<s.det());return tmp1;}
int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
if(fmin(rn,rm)<fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
tmp1(0,0)=0;
return tmp1;
}
template<class T>MATRIX<T> operator>(MATRIX<T> r,MATRIX<T> s)
{
MATRIX<T> tmp1(1);
if(r.size()>s.size())
{tmp1(0,0)=1;return tmp1;}
else if(r.size()<s.size())
{tmp1(0,0)=0;return tmp1;}
else if((r.size_M().n==r.size_M().m)&(s.size_M().n==r.size_M().m))
{tmp1(0,0)=(T)(r.det()>s.det());return tmp1;}
int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
if(fmin(rn,rm)>fmin(sn,sm)){tmp1(0,0)=1;return tmp1;}
tmp1(0,0)=0;
return tmp1;
}
//template<class T>MATRIX<T> operator==(MATRIX<T> r,MATRIX<T> s) //définie par défaut
// {
// MATRIX<T> tmp1(1);
// int rn=r.size_M().n,rm=r.size_M().m,sn=s.size_M().n,sm=s.size_M().m;
// if((rn!=sn)|(rm!=sm)){tmp1(0,0)=0;return tmp1;}
// else
// for(int i=0;i<rn;i++)
// for(int j=0;j<rm;j++)
// if(r(i,j)!=s(i,j)){tmp1(0,0)=0;return tmp1;}
// tmp(0,0)=1;
// return tmp1;
// }
#endif //ndef _matrice
/*int diagonalisation(MATRIX2,MATRIX2&,MATRIX2&,std::vector<double>&);
int diagonalisation(MATRIX2 A,MATRIX2& D,MATRIX2& P,std::vector<double>& Lambda)
{
double s,df;
D=MATRIX2(A.n);P=MATRIX2(A.n);Lambda.clear();
if (A.n!=A.m) {__ENVOI("imposible de diagonalise une matrice non carre!\n");exit(0);}
double inter=0,tmp1,tmp2,tmp3;
double pr=200;//precision
for (int i=0;i<A.n*A.n;i++)inter+=fabs(A.M[i]);
for (double l=-inter;l<inter;l+=inter/pr)
/ *{
while(1)
{
df=((A-l*MATRIX2(A.n).Id()).det()-(A-(l-inter/(20*pr))*MATRIX2(A.n).Id()).det())/(inter/(20*pr));
s=l+(A-l*MATRIX2(A.n).Id()).det()/df;cout<<s-l<<endl;if(s<l) {l+=inter;continue;}
if (s-l!=0) l=s;
else break;if (!-inter<l) system("pause");
}
Lambda.push_back(l);
}//méthode de Newton(proto necessite d'avoir la dérivée du polynome caractéristique)* /
/ * {
tmp1=(A-l*MATRIX2(A.n).Id()).det();
if (tmp1*tmp2<0 & l!=-inter)
Lambda.push_back(l-inter*2.0/pr);
if (tmp1==0)
Lambda.push_back(l);
tmp2=tmp1;
} //méthode perso, lente, possibilité de sauter des solutions, et impresice
if (Lambda.size()!=A.n) {cout<<MATRIX2(Lambda);__ENVOI("probleme de nombre de VP\n");exit(0);}
for (int i=0;i<A.n;i++)D(i,i)=Lambda[i];
tmp2=0;for (int i=1;i<A.n;i++) {if (A(0,i)!=0)tmp2=1;}tmp3=tmp2;//traitement du cas où la 1ére coordonnée du VP est nulle.(proto)
MATRIX2 AV(A.n-1),AVt;
for(int i=0;i<A.n;i++)
{
for (int k=0;k<A.n-1;k++)for (int j=0;j<A.n-1;j++)AV(k,j)=A(k+1,j+1);
AV=AV-Lambda[i]*MATRIX2(AV.n).Id();
if(abs(A(0,0)-Lambda[i])<2*inter/pr) tmp3=1;else tmp3=tmp2;
tmp1=AV.det();
P(0,i)=tmp3;
for(int j=1;j<A.n;j++)
{
AVt=AV;
for (int k=0;k<AV.n;k++) AVt(k,j-1)=-A(k+1,0)/ **tmp3* /;
/ * P(j,i)=AVt.det()/tmp1;
}
}
return 0;
}* /
/ * pour passer n m et M en private, il faut definir les fonction qui les utilise en friend dans la declaration
(si j'ai bien compris). et rajouter eventuellement quelques fonctions, comme resize pour aller avec la fonction
push_back, si l'on ne veut pas une matrice carré */
Sources de la même categorie
Commentaires et avis
Discussions en rapport avec ce code source dans le forum
Template matrice [ par anisdilou ]
Je souhaite q'on s'aide a fin d'meliorer nos nivauxSalut : Je veux implémenter une classe qui représente des matrices de dimensions et type de donnés
Ouverture d'une image .bmp en tant que matrice [ par jpout ]
Bonjour,Je cherche à ouvrir une image (*bmp ou *.jpg) en tant que simple matrice afin de pouvoir travailler plus facilement sur les pixels. Je cherche
opérations matrice /temps exécution prg [ par 0wil0 ]
Bonjour, J'effectue dans mon programme des opérations relativement simples sur des matrices (additions, soustractions, moyenne des éléments de matrice
valeurs propres et vecteurs associés [ par GUARMAH ]
Bonjour,je cherche un code source qui prend en entrée une matrice réelle plein qcq et qui donne TOUTES les valeurs propres ainsi que les vecteurs prop
Pb quaternion [ par Arnaud16022 ]
alors voila, la pb est de faire tourner un point autour d'une ligne.comme je ne suis pas le premier a ma demander comment faire, j'ai cherché un peu e
MULTIPLICATION D'UNE MATRICE PAR UN VECTEUR [ par jfk20004 ]
Quelqu'un pourrait il m'expliquer le bout de code suivant tiré d'un prog de raytrace .Cette partie est censée multiplier une matriceet un vecteur.Je n
Appliquer une matrice 3*3 à image en C++ [ par dj roupe ]
Je dispose d'une image définie par un ensemble de pixels et je souhaite appliquer une matrice homogène 3*3 à cette image et je ne sais pas comment fai
affichage 3d [ par Arnaud16022 ]
voici le probleme:Soit A un point défini par la ctruct Vecteur3d{float x,x,z;} de coordonées dans l'espace cartésien orienté orthonormé (A.x,A.y,A.z).
Multiplication d'une matrice par un scalaire [ par skrime ]
Bonjour, j'ai un exercice à faire qui consiste à multiplier une matrice par un scalaire en C (la prof ne veut pas qu'on se serve des boucles FOR), je
comment sauvegarder une matrice sous fichier et la racharger qu'on veut?? [ par malbb2000 ]
salut tt le monde j'aimerais bien que qqn puisse m'aider sur le sujet que j'ai posais d'avance merci en core
|
Derniers Blogs
COMMENT MAPPER UNE VUE SQL SUR UNE COLLECTION DE COMPLEX TYPE?COMMENT MAPPER UNE VUE SQL SUR UNE COLLECTION DE COMPLEX TYPE? par Matthieu MEZIL
Avec EF, les vues doivent être mappées sur des entity types. Le problème c'est que les entity types doivent avoir une clé. Avec EF, nous avons les complex type qui n'ont pas de clé mais les vues ne peuvent pas être mappées dessus. Avec EF4, il est possibl...
Cliquez pour lire la suite de l'article par Matthieu MEZIL [WF4] UN BINDING ACTIVITY/ACTIVITYDESIGNER QUI PASSE MAL?[WF4] UN BINDING ACTIVITY/ACTIVITYDESIGNER QUI PASSE MAL? par JeremyJeanson
Certain d'entre vous on peut être vécu cette situation embarrassante après quelques temps passer avec WF4 : Au début avec mon " ActivityDesigner" , tout allait bien. Et puis un jour j'ai au des problèmes de " Binding" . Alors nous sommes allé sur le site ...
Cliquez pour lire la suite de l'article par JeremyJeanson MYTIC - SHAREPOINT 2010 : DéJà UN MYTHE MICROSOFT ?MYTIC - SHAREPOINT 2010 : DéJà UN MYTHE MICROSOFT ? par junarnoalg
La prochaine session de MyTIC aura lieu à Namur, le 23 mars prochain. Pendant presque une heure, nous parlerons de SharePoint 2010. Voici un aperçu du programme.
Accueil : 17h30 Début de la session : 18h00 - Les nouvelles int...
Cliquez pour lire la suite de l'article par junarnoalg
Forum
ERREUR DE POINTEURERREUR DE POINTEUR par africanwinners
Cliquez pour lire la suite par africanwinners CLISTCTRLCLISTCTRL par dorras7
Cliquez pour lire la suite par dorras7
Logiciels
Academy System (10.9.4.0)ACADEMY SYSTEM (10.9.4.0)Logiciel de gestion des établissements.
- élèves/étudiants (inscription, dossier, absence...)
-... Cliquez pour télécharger Academy System Xilisoft Convertisseur Vidéo Ultimate (5.1.39.0305)XILISOFT CONVERTISSEUR VIDéO ULTIMATE (5.1.39.0305)Xilisoft Convertisseur Vidéo Ultimate est un outil puissant de conversion vidéo, facile à utilise... Cliquez pour télécharger Xilisoft Convertisseur Vidéo Ultimate Xilisoft DVD Ripper Ultimate (5.0.64.0304)XILISOFT DVD RIPPER ULTIMATE (5.0.64.0304)Xilisoft DVD Ripper Ultimate est un logiciel excellent pour copier et convertir DVD vers presque ... Cliquez pour télécharger Xilisoft DVD Ripper Ultimate Rigs of Rods (63.3)RIGS OF RODS (63.3)c'est un jeu de multi-simulation camions,autobus voitures, avions, bateaux, hélicoptère avec défo... Cliquez pour télécharger Rigs of Rods
|