Class matrice avec template

Soyez le premier à donner votre avis sur cette source.

Vue 7 528 fois - Téléchargée 301 fois

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 / Exemple :


/* 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é */

Codes Sources

A voir également

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.