Index: /trunk/SophyaLib/TArray/basarr.cc
===================================================================
--- /trunk/SophyaLib/TArray/basarr.cc	(revision 2582)
+++ /trunk/SophyaLib/TArray/basarr.cc	(revision 2583)
@@ -29,6 +29,48 @@
 short  BaseArray::default_vector_type = ColumnVector;
 sa_size_t BaseArray::openmp_size_threshold = 200000;
+uint_2 BaseArray::matrix_product_optim = 1;
 
 // ------ Methodes statiques globales --------
+
+//! Set optimization strategy for matrix product
+/*!
+  \param opt : bit coded optimization strategy
+  \warning Default is opt=1
+  \verbatim
+  bit 0 : choose matrix product optimization or not
+        0=no optimization, 1=optimization
+  bit 1 : force optimization in any case (only if bit0=1)
+          (if not the TMatrix::Multiply method decide or
+           not if the product should be optimized)
+        0=do not force and let the method decide, 1=force
+  bit 2 : optimize the product A * B when A-columns and
+          B-rows are not packed (for example if the product
+          is C = ("A Fortran-like" * "B C-like").
+        . That will do a copy of one of the two matrices,
+          so that will result in an increase of the memory
+          space needed.
+        . For big matrices that decrease the computing time.
+        . Do not use this optimisation for small matrices
+          because that would increase the computing time.
+        0=do not optimze that way, 1=optimze that way
+  \endverbatim
+  \verbatim
+   Sumary of the allowed values for "opt"
+           0 = no optimization at all (whatever the other bits are)
+           1 = optimize but let the method decides if optimization
+               is needed regarding the sizes of the matrices.
+           3 = force optimization whatever the sizes of the matrices are.
+           5 = optimisation with method decision ("1")
+               AND optimize by copying when A-columns and B-rows
+               are not packed
+           7 = force optimization ("3")
+               AND optimize by copying when A-columns and B-rows
+               are not packed
+  \endverbatim
+*/
+void BaseArray::SetMatProdOpt(uint_2 opt)
+{
+  matrix_product_optim = opt;
+}
 
 //! Set maximum number of printed elements and print level
Index: /trunk/SophyaLib/TArray/basarr.h
===================================================================
--- /trunk/SophyaLib/TArray/basarr.h	(revision 2582)
+++ /trunk/SophyaLib/TArray/basarr.h	(revision 2583)
@@ -56,4 +56,8 @@
   //! Get Default Vector Type
   static inline short GetDefaultVectorType() { return default_vector_type; }
+
+  //! Optimization choice for matrix product
+  static void SetMatProdOpt(uint_2 opt=1);
+  static inline uint_2 GetMatProdOpt(void) {return matrix_product_optim;}
   
   // Creator / destructor
@@ -199,4 +203,5 @@
   static short default_vector_type;    //!< Default vector type Row/Column
   static sa_size_t openmp_size_threshold; //!< Size limit for parallel routine calls
+  static uint_2 matrix_product_optim; //!< optimization level for matrix product
 };
 
Index: /trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.cc	(revision 2582)
+++ /trunk/SophyaLib/TArray/tmatrix.cc	(revision 2583)
@@ -1,3 +1,3 @@
-// $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 +6,4 @@
 #include "pexceptions.h"
 #include "tmatrix.h"
-
-//#define DO_NOT_OPTIMIZE_PRODUCT
 
 /*!
@@ -448,12 +446,14 @@
 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,21 +461,34 @@
   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,7 +506,13 @@
       }
     }
-  }
-  // 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,7 +529,13 @@
     }
     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,7 +552,13 @@
     }
     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,7 +575,35 @@
     }
     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,7 +620,10 @@
     }
     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,10 +640,13 @@
     }
     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;
 }
