Changeset 2583 in Sophya


Ignore:
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

Location:
trunk/SophyaLib/TArray
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/basarr.cc

    r2412 r2583  
    2929short  BaseArray::default_vector_type = ColumnVector;
    3030sa_size_t BaseArray::openmp_size_threshold = 200000;
     31uint_2 BaseArray::matrix_product_optim = 1;
    3132
    3233// ------ Methodes statiques globales --------
     34
     35//! Set optimization strategy for matrix product
     36/*!
     37  \param opt : bit coded optimization strategy
     38  \warning Default is opt=1
     39  \verbatim
     40  bit 0 : choose matrix product optimization or not
     41        0=no optimization, 1=optimization
     42  bit 1 : force optimization in any case (only if bit0=1)
     43          (if not the TMatrix::Multiply method decide or
     44           not if the product should be optimized)
     45        0=do not force and let the method decide, 1=force
     46  bit 2 : optimize the product A * B when A-columns and
     47          B-rows are not packed (for example if the product
     48          is C = ("A Fortran-like" * "B C-like").
     49        . That will do a copy of one of the two matrices,
     50          so that will result in an increase of the memory
     51          space needed.
     52        . For big matrices that decrease the computing time.
     53        . Do not use this optimisation for small matrices
     54          because that would increase the computing time.
     55        0=do not optimze that way, 1=optimze that way
     56  \endverbatim
     57  \verbatim
     58   Sumary of the allowed values for "opt"
     59           0 = no optimization at all (whatever the other bits are)
     60           1 = optimize but let the method decides if optimization
     61               is needed regarding the sizes of the matrices.
     62           3 = force optimization whatever the sizes of the matrices are.
     63           5 = optimisation with method decision ("1")
     64               AND optimize by copying when A-columns and B-rows
     65               are not packed
     66           7 = force optimization ("3")
     67               AND optimize by copying when A-columns and B-rows
     68               are not packed
     69  \endverbatim
     70*/
     71void BaseArray::SetMatProdOpt(uint_2 opt)
     72{
     73  matrix_product_optim = opt;
     74}
    3375
    3476//! Set maximum number of printed elements and print level

  • trunk/SophyaLib/TArray/basarr.h

    r2322 r2583  
    5656  //! Get Default Vector Type
    5757  static inline short GetDefaultVectorType() { return default_vector_type; }
     58
     59  //! Optimization choice for matrix product
     60  static void SetMatProdOpt(uint_2 opt=1);
     61  static inline uint_2 GetMatProdOpt(void) {return matrix_product_optim;}
    5862 
    5963  // Creator / destructor
    … …  
    199203  static short default_vector_type;    //!< Default vector type Row/Column
    200204  static sa_size_t openmp_size_threshold; //!< Size limit for parallel routine calls
     205  static uint_2 matrix_product_optim; //!< optimization level for matrix product
    201206};
    202207

  • trunk/SophyaLib/TArray/tmatrix.cc

    r2575 r2583  
    1 // $Id: tmatrix.cc,v 1.26 2004-07-29 12:31:16 ansari Exp $
     1// $Id: tmatrix.cc,v 1.27 2004-07-30 10:24:12 cmv Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    … …  
    66#include "pexceptions.h"
    77#include "tmatrix.h"
    8 
    9 //#define DO_NOT_OPTIMIZE_PRODUCT
    108
    119/*!
    … …  
    448446template <class T>
    449447TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
     448// Calcul de C= rm = A*B   (A=*this)
     449// Remember: C-like      matrices are column packed
     450//           Fortan-like matrices are line   packed
    450451{
    451452  if (NCols() != b.NRows())
    452453    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
    453454
    454   // Commentaire: pas de difference de vitesse notable selon les mappings choisis pour la matrice produit "rm"
     455  // Commentaire: pas de difference de vitesse notable selon le mapping de la matrice produit "rm"
    455456  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
    456457  TMatrix<T> rm(NRows(), b.NCols(), mm);
    457   rm = (T) 0;  // Rien ne garantit que rm soit mise a zero a la creation !
    458458
    459459  // Les "steps" pour l'adressage des colonnes de A et des lignes de B
    … …  
    461461  sa_size_t stepb = b.Step(b.RowsKA());
    462462
    463   // On decide si on optimise ou non selon les dimensions de A et B
    464   // (il semble que optimiser ou non ne degrade pas notablement la vitesse pour les petites matrices)
    465 #ifdef DO_NOT_OPTIMIZE_PRODUCT
    466   bool no_optim = true;
    467 #else
    468   bool no_optim = false;
    469   //NON juste pour memoire: if((stepa+stepb)*NCols()*sizeof(T)<100000) no_optim=true;
    470 #endif
    471 
    472   // Calcul de C=rm = A*B   (A=*this)
    473   // Remember: C-like      matrices are column packed
    474   //           Fortan-like matrices are line   packed
     463  // Taille totale des matrices A et B
     464  size_t totsiza = this->DataBlock().Size();
     465  size_t totsizb = b.DataBlock().Size();
     466
     467
     468  ///////////////////////////////////////////////////////////////////
     469  // On decide si on optimise ou non selon les dimensions de A et B //
     470  // (il semble que optimiser ou non ne degrade pas                 //
     471  //  beaucoup la vitesse pour les petites matrices)                //
     472  ////////////////////////////////////////////////////////////////////
     473
     474  uint_2 popt = GetMatProdOpt();
     475  bool no_optim = false; // optimization demandee par default
     476  if( (popt&(uint_2)1) == 0 ) { // pas d'optimization explicitement demande
     477    no_optim = true;
     478  } else if( (popt&(uint_2)2) == 0 ) { // pas d'optimization forcee, la methode decide
     479    // On part sur une disponibilite dans le cache processeur de 100 ko
     480    // (A et B peuvent etre stoquees dans le cache)
     481    if((totsiza+totsizb)*sizeof(T)<100000) no_optim = true;
     482  }
     483
    475484  sa_size_t r,c,k;
    476485  T sum;
    477486  const T * pe;
    478487
    479   // Pas d'optimisation demandee
     488
     489  /////////////////////////////////
     490  // Pas d'optimisation demandee //
     491  /////////////////////////////////
     492
    480493  if( no_optim ) {
    481494    //cout<<"no_optim("<<no_optim<<") "<<stepa<<" "<<stepb<<endl;
    … …  
    493506      }
    494507    }
    495   }
    496   // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme)
    497   else if(stepa==1 && stepb==1) {
     508    return rm;
     509  }
     510
     511
     512  //////////////////////////////////////////////////////////////////////////////////
     513  // A.col est packed  et  B.row est packed (on a interet a optimiser quand meme) //
     514  //////////////////////////////////////////////////////////////////////////////////
     515
     516  if(stepa==1 && stepb==1) {
    498517    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
    499518     T * pea = new T[rm.NCols()];
    … …  
    510529    }
    511530    delete [] pea;
    512  }
    513   // A.col est packed  et  B.row n'est pas packed
    514   else if(stepa==1 && stepb!=1) {
     531    return rm;
     532  }
     533
     534
     535  //////////////////////////////////////////////////
     536  // A.col est packed  et  B.row n'est pas packed //
     537  //////////////////////////////////////////////////
     538
     539  if(stepa==1 && stepb!=1) {
    515540    //cout<<"A.col packed && B.row not packed "<<stepa<<" "<<stepb<<endl;
    516541    const T * pea;
    … …  
    527552    }
    528553    delete [] peb;
    529   }
    530   // A.col n'est pas packed  et  B.row est packed
    531   else if(stepa!=1 && stepb==1) {
     554    return rm;
     555  }
     556
     557
     558  //////////////////////////////////////////////////
     559  // A.col n'est pas packed  et  B.row est packed //
     560  //////////////////////////////////////////////////
     561
     562  if(stepa!=1 && stepb==1) {
    532563    //cout<<"A.col not packed && B.row packed "<<stepa<<" "<<stepb<<endl;
    533564    T * pea = new T[rm.NCols()];
    … …  
    544575    }
    545576    delete [] pea;
    546   }
    547   // A.col n'est pas packed  et  B.row n'est pas packed, stepb>stepa
    548   else if(stepa!=1 && stepb!=1 && stepb>stepa) {
     577    return rm;
     578  }
     579
     580
     581  ////////////////////////////////////////////////////////
     582  // A.col n'est pas packed  et  B.row n'est pas packed //
     583  ////////////////////////////////////////////////////////
     584
     585  //---- On demande l'optimization par copie d'une des matrices
     586
     587  if( (popt&(uint_2)4) ) {
     588    // On copie la plus petite
     589    if(totsiza<totsizb) {   // on copie A
     590      //cout<<"A.col not packed && B.row not packed ==> copy A to optimize "<<stepa<<" "<<stepb<<endl;
     591      // Acopy doit etre C-like pour etre column-packed
     592      TMatrix<T> acopy(NRows(),NCols(),BaseArray::CMemoryMapping);
     593      acopy = *this;
     594      rm = acopy.Multiply(b,mm);
     595    } else {               // on copie B
     596      //cout<<"A.col not packed && B.row not packed ==> copy B to optimize "<<stepa<<" "<<stepb<<endl;
     597      // Bcopy doit etre Fortran-like pour etre column-packed
     598      TMatrix<T> bcopy(b.NRows(),b.NCols(),BaseArray::FortranMemoryMapping);
     599      bcopy = b;
     600      rm = Multiply(bcopy,mm);
     601    }
     602    return rm;
     603  }
     604 
     605  //---- stepb>stepa
     606
     607  if(stepa!=1 && stepb!=1 && stepb>stepa) {
    549608    //cout<<"A.col not packed && B.row not packed ==> optimize on B "<<stepa<<" "<<stepb<<endl;
    550609    const T * pea;
    … …  
    561620    }
    562621    delete [] peb;
    563   }
    564   // A.col n'est pas packed  et  B.row n'est pas packed, stepa>=stepb
    565   else if(stepa!=1 && stepb!=1) {
     622    return rm;
     623  }
     624
     625  //---- stepa>=stepb
     626
     627  if(stepa!=1 && stepb!=1) {
    566628    //cout<<"A.col not packed && B.row not packed ==> optimize on A "<<stepa<<" "<<stepb<<endl;
    567629    T * pea = new T[rm.NCols()];
    … …  
    578640    }
    579641    delete [] pea;
    580   }
    581   else {
    582     cout<<"TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! "<<endl;
    583     throw(SzMismatchError("TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! ") );
    584   }
    585 
     642    return rm;
     643  }
     644
     645
     646  //////////////////////////////////////////////////
     647  // Cas non prevu, on ne doit JAMAIS arriver ici //
     648  //////////////////////////////////////////////////
     649  cout<<"TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! "<<endl;
     650  throw(SzMismatchError("TMatrix<T>::Multiply(b) Optimize case not treated... Please report BUG !!! ") );
    586651  return rm;
    587652}
Note: See TracChangeset for help on using the changeset viewer.