Changeset 2583 in Sophya for trunk/SophyaLib/TArray/tmatrix.cc

Timestamp:

Jul 30, 2004, 12:24:12 PM (21 years ago)

Author:

cmv

Message:

• Intro decision auto d'optimisation produit de matrices
• Possibilite a l'utilisateur pour choisir l'optimisation
• cas FxC optimise par copie

cmv 30/07/04

File:

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

trunk/SophyaLib/TArray/tmatrix.cc

@@ -1 @@
-// $Id: tmatrix.cc,v 1.26 2004-07-29 12:31:16 ansari Exp $
+// $Id: tmatrix.cc,v 1.27 2004-07-30 10:24:12 cmv Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -6 +6 @@
 #include "pexceptions.h"
 #include "tmatrix.h"
-
-//#define DO_NOT_OPTIMIZE_PRODUCT
 
 /*!
@@ -448 +446 @@
 template <class T>
 TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
+// Calcul de C= rm = A*B   (A=*this)
+// Remember: C-like      matrices are column packed
+//           Fortan-like matrices are line   packed
 {
   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"
+  // Commentaire: pas de difference de vitesse notable selon le mapping de 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
@@ -461 +461 @@
   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
+  // Taille totale des matrices A et B
+  size_t totsiza = this->DataBlock().Size();
+  size_t totsizb = b.DataBlock().Size();
+
+
+  ///////////////////////////////////////////////////////////////////
+  // On decide si on optimise ou non selon les dimensions de A et B //
+  // (il semble que optimiser ou non ne degrade pas                 //
+  //  beaucoup la vitesse pour les petites matrices)                //
+  ////////////////////////////////////////////////////////////////////
+
+  uint_2 popt = GetMatProdOpt();
+  bool no_optim = false; // optimization demandee par default
+  if( (popt&(uint_2)1) == 0 ) { // pas d'optimization explicitement demande
+    no_optim = true;
+  } else if( (popt&(uint_2)2) == 0 ) { // pas d'optimization forcee, la methode decide
+    // On part sur une disponibilite dans le cache processeur de 100 ko
+    // (A et B peuvent etre stoquees dans le cache)
+    if((totsiza+totsizb)*sizeof(T)<100000) no_optim = true;
+  }
+
   sa_size_t r,c,k;
   T sum;
   const T * pe;
 
-  // Pas d'optimisation demandee
+
+  /////////////////////////////////
+  // Pas d'optimisation demandee //
+  /////////////////////////////////
+
   if( no_optim ) {
     //cout<<"no_optim("<<no_optim<<") "<<stepa<<" "<<stepb<<endl;
@@ -493 +506 @@
       }
     }
-  }
-  // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme)
-  else if(stepa==1 && stepb==1) {
+    return rm;
+  }
+
+
+  //////////////////////////////////////////////////////////////////////////////////
+  // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme) //
+  //////////////////////////////////////////////////////////////////////////////////
+
+  if(stepa==1 && stepb==1) {
     //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
      T * pea = new T[rm.NCols()];
@@ -510 +529 @@
     }
     delete [] pea;
- }
-  // A.col est packed  et  B.row n'est pas packed
-  else if(stepa==1 && stepb!=1) {
+    return rm;
+  }
+
+
+  //////////////////////////////////////////////////
+  // A.col est packed  et  B.row n'est pas packed //
+  //////////////////////////////////////////////////
+
+  if(stepa==1 && stepb!=1) {
     //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
     const T * pea;
@@ -527 +552 @@
     }
     delete [] peb;
-  }
-  // A.col n'est pas packed  et  B.row est packed
-  else if(stepa!=1 && stepb==1) {
+    return rm;
+  }
+
+
+  //////////////////////////////////////////////////
+  // A.col n'est pas packed  et  B.row est packed //
+  //////////////////////////////////////////////////
+
+  if(stepa!=1 && stepb==1) {
     //cout<<"A.col not packed && B.row packed "<<stepa<<" "<<stepb<<endl;
     T * pea = new T[rm.NCols()];
@@ -544 +575 @@
     }
     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) {
+    return rm;
+  }
+
+
+  ////////////////////////////////////////////////////////
+  // A.col n'est pas packed  et  B.row n'est pas packed //
+  ////////////////////////////////////////////////////////
+
+  //---- On demande l'optimization par copie d'une des matrices
+
+  if( (popt&(uint_2)4) ) {
+    // On copie la plus petite
+    if(totsiza<totsizb) {   // on copie A
+      //cout<<"A.col not packed && B.row not packed ==> copy A to optimize "<<stepa<<" "<<stepb<<endl;
+      // Acopy doit etre C-like pour etre column-packed
+      TMatrix<T> acopy(NRows(),NCols(),BaseArray::CMemoryMapping);
+      acopy = *this;
+      rm = acopy.Multiply(b,mm);
+    } else {               // on copie B
+      //cout<<"A.col not packed && B.row not packed ==> copy B to optimize "<<stepa<<" "<<stepb<<endl;
+      // Bcopy doit etre Fortran-like pour etre column-packed
+      TMatrix<T> bcopy(b.NRows(),b.NCols(),BaseArray::FortranMemoryMapping);
+      bcopy = b;
+      rm = Multiply(bcopy,mm);
+    }
+    return rm;
+  }
+ 
+  //---- stepb>stepa
+
+  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;
@@ -561 +620 @@
     }
     delete [] peb;
-  }
-  // A.col n'est pas packed  et  B.row n'est pas packed, stepa>=stepb
-  else if(stepa!=1 && stepb!=1) {
+    return rm;
+  }
+
+  //---- stepa>=stepb
+
+  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()];
@@ -578 +640 @@
     }
     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;
+  }
+
+
+  //////////////////////////////////////////////////
+  // Cas non prevu, on ne doit JAMAIS arriver ici //
+  //////////////////////////////////////////////////
+  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;
 }