Changeset 2574 in Sophya for trunk/SophyaLib/TArray

Timestamp:

Jul 29, 2004, 10:40:49 AM (21 years ago)

Author:

cmv

Message:

Optimisation produit de matrice a la sauce cmv 29/07/04

File:

• trunk/SophyaLib/TArray/tmatrix.cc (modified) (4 diffs)

trunk/SophyaLib/TArray/tmatrix.cc

@@ -1 @@
-// $Id: tmatrix.cc,v 1.24 2003-08-08 13:07:06 ansari Exp $
+// $Id: tmatrix.cc,v 1.25 2004-07-29 08:40:49 cmv Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -6 +6 @@
 #include "pexceptions.h"
 #include "tmatrix.h"
+
+//#define DO_NOT_OPTIMIZE_PRODUCT
 
 /*!
@@ -411 +413 @@
   \param mm : define the memory mapping type for the return matrix
  */
+////////////// Routine de base sans optimisation //////////////
+/*
 template <class T>
 TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
@@ -437 +441 @@
   return rm;
 }
+*/
+
+////////////// Routine optimisee //////////////
+template <class T>
+TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
+{
+  if (NCols() != b.NRows()) 
+    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
+
+  // Commentaire: pas de difference de vitesse notable selon les mappings choisis pour la matrice produit "rm"
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> rm(NRows(), b.NCols(), mm);
+  rm = (T) 0;  // Rien ne garantit que rm soit mise a zero a la creation !
+
+  // Les "steps" pour l'adressage des colonnes de A et des lignes de B
+  sa_size_t stepa = Step(ColsKA());
+  sa_size_t stepb = b.Step(b.RowsKA());
+
+  // On decide si on optimise ou non selon les dimensions de A et B
+  // (il semble que optimiser ou non ne degrade pas notablement la vitesse pour les petites matrices)
+#ifdef DO_NOT_OPTIMIZE_PRODUCT
+  bool no_optim = true;
+#else
+  bool no_optim = false;
+  //NON juste pour memoire: if((stepa+stepb)*NCols()*sizeof(T)<100000) no_optim=true;
+#endif
+
+  // Calcul de C=rm = A*B   (A=*this)
+  // Remember: C-like      matrices are column packed
+  //           Fortan-like matrices are line   packed
+  sa_size_t r,c,k;
+  T sum;
+  const T * pe;
+
+  // Pas d'optimisation demandee
+  if( no_optim ) {
+    //cout<<"no_optim("<<no_optim<<") "<<stepa<<" "<<stepb<<endl;
+    const T * pea;
+    const T * peb;
+    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
+      for(c=0; c<rm.NCols(); c++) {   // Boucle sur les colonnes de B
+        sum = 0;
+        pea = &((*this)(r,0));
+        peb = &(b(0,c));
+        // On gagne un peu en remplacant "pea[k*stepa]" par "pea+=stepa" pour les grosses matrices
+        //for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
+        for(k=0; k<NCols(); k++) {sum += (*pea)*(*peb); pea+=stepa; peb+=stepb;}
+        rm(r,c) = sum;
+      }
+    }
+  }
+  // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme)
+  else if(stepa==1 && stepb==1) {
+    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
+     T * pea = new T[rm.NCols()];
+    const T * peb;
+    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
+      pe = &((*this)(r,0));
+      for(k=0; k<NCols(); k++) {pea[k] = *(pe++);}
+      for(c=0; c<rm.NCols(); c++) {   // Boucle sur les colonnes de B
+        sum = 0;
+        peb = &(b(0,c));
+        for(k=0; k<NCols(); k++) sum += *(peb++)*pea[k];
+        rm(r,c) = sum;
+      }
+    }
+    delete [] pea;
+ }
+  // A.col est packed  et  B.row n'est pas packed
+  else if(stepa==1 && stepb!=1) {
+    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
+    const T * pea;
+    T * peb = new T[rm.NCols()];
+    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
+      pe = &(b(0,c));
+      for(k=0; k<NCols(); k++) {peb[k] = *pe; pe+=stepb;}
+      for(r=0; r<rm.NRows(); r++) {   // Boucle sur les lignes de A
+        sum = 0;
+        pea = &((*this)(r,0));
+        for(k=0; k<NCols(); k++) sum += pea[k]*peb[k];
+        rm(r,c) = sum;
+      }
+    }
+    delete [] peb;
+  }
+  // A.col n'est pas packed  et  B.row est packed
+  else if(stepa!=1 && stepb==1) {
+    //cout<<"A.col not packed && B.row packed "<<stepa<<" "<<stepb<<endl;
+    T * pea = new T[rm.NCols()];
+    const T * peb;
+    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
+      pe = &((*this)(r,0));
+      for(k=0; k<NCols(); k++) {pea[k] = *pe; pe+=stepa;}
+      for(c=0; c<rm.NCols(); c++) {   // Boucle sur les colonnes de B
+        sum = 0;
+        peb = &(b(0,c));
+        for(k=0; k<NCols(); k++) sum += pea[k]*peb[k];
+        rm(r,c) = sum;
+      }
+    }
+    delete [] pea;
+  }
+  // A.col n'est pas packed  et  B.row n'est pas packed, stepb>stepa
+  else if(stepa!=1 && stepb!=1 && stepb>stepa) {
+    //cout<<"A.col not packed && B.row not packed ==> optimize on B "<<stepa<<" "<<stepb<<endl;
+    const T * pea;
+    T * peb = new T[rm.NCols()];
+    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
+      pe = &(b(0,c));
+      for(k=0; k<NCols(); k++) {peb[k] = *pe; pe+=stepb;}
+      for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
+        sum = 0;
+        pea = &((*this)(r,0));
+        for(k=0; k<NCols(); k++) {sum += (*pea)*peb[k]; pea+=stepa;}
+        rm(r,c) = sum;
+      }
+    }
+    delete [] peb;
+  }
+  // A.col n'est pas packed  et  B.row n'est pas packed, stepa>=stepb
+  else if(stepa!=1 && stepb!=1) {
+    //cout<<"A.col not packed && B.row not packed ==> optimize on A "<<stepa<<" "<<stepb<<endl;
+    T * pea = new T[rm.NCols()];
+    const T * peb;
+    for(r=0; r<rm.NRows(); r++) {     // Boucle sur les lignes de A
+      pe = &((*this)(r,0));
+      for(k=0; k<NCols(); k++) {pea[k] = *pe; pe+=stepa;}
+      for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
+        sum = 0;
+        peb = &(b(0,c));
+        for(k=0; k<NCols(); k++) {sum += pea[k]*(*peb); peb+=stepb;}
+        rm(r,c) = sum;
+      }
+    }
+    delete [] pea;
+  }
+  else {
+    cout<<"TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! "<<endl;
+    throw(SzMismatchError("TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! ") );
+  }
+
+  return rm;
+}
 
 ///////////////////////////////////////////////////////////////