Changeset 804 in Sophya


Ignore:
Timestamp:
Apr 3, 2000, 7:36:01 PM (25 years ago)
Author:
ansari
Message:

Amelioation / debugging de la classe TArray<T> - TVector et TMatrix

heritent maintenant de TArray<T> - Classe RCMatrix rendu prive au fichier
sopemtx.cc - linfit.cc integre a sopemtx.cc

Reza 03/04/2000

Location:
trunk/SophyaLib/TArray
Files:
1 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/basarr.cc

    r787 r804  
    1414                                 "Stride(int )", "ElemCheckBound()", "operator()" };
    1515uint_4 BaseArray::max_nprt_ = 50;
     16short  BaseArray::default_memory_mapping = CMemoryMapping;
    1617
    1718// Methodes statiques globales
    … …  
    2223
    2324
    24 uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset)
     25uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset)
    2526{
    2627  uint_8 rs = step;
    … …  
    2829  return(rs+offset);
    2930}
     31
     32short BaseArray::SetDefaultMemoryMapping(short mm)
     33{
     34  default_memory_mapping = ( (mm == CMemoryMapping) ? CMemoryMapping : FortranMemoryMapping) ;
     35  return default_memory_mapping;
     36}
     37
     38
     39short BaseArray::SelectMemoryMapping(short mm)
     40// si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
     41// sinon (mm == FortranMemoryMapping)  elements d'une colonne suivant X
     42{
     43  if ( (mm == CMemoryMapping) || (mm == FortranMemoryMapping) )  return (mm) ;
     44  else return (default_memory_mapping);
     45}
     46
     47void BaseArray::UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny)
     48{
     49  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
     50  if (mm == CMemoryMapping) {
     51    marowi_ = veceli_ = 1;  macoli_ = 0;
     52    uint_4 nr = nx;
     53    nx = ny;   ny = nr;
     54  }
     55  else {
     56    marowi_ = veceli_ = 0;  macoli_ = 1;
     57  } 
     58}
     59
     60void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
     61{
     62  if (mm == SameMemoryMapping)
     63    mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
     64  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
     65  if (mm == CMemoryMapping) {
     66    marowi_ = veceli_ = 1;  macoli_ = 0;
     67  }
     68  else {
     69    marowi_ = veceli_ = 0;  macoli_ = 1;
     70  } 
     71}
     72
    3073
    3174// -------------------------------------------------------
    … …  
    4487  moystep_ = 0;
    4588  offset_ = 0;
    46   // Default for matrices : Lines = Y , Columns = X
     89  // Default for matrices : Lines = Y , Columns = X, Default Vectore= ColumnVector
    4790  macoli_ = 0;
    48   marowi_ = 1;
     91  veceli_ = marowi_ = 1;
    4992}
    5093
    … …  
    60103  for(int k=0; k<ndim_; k++)
    61104    if (size_[k] != a.size_[k])  return(false);
    62   if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) )  return(false);
     105  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
     106       (veceli_ != a.veceli_) )  return(false);
    63107  return(true);
    64108}
    … …  
    196240  // Default for matrices : Lines = Y , Columns = X
    197241  macoli_ = 0;
    198   marowi_ = 1;
     242  marowi_ = veceli_ = 1;
    199243  // Update OK
    200244  ndim_ = ndim;
    … …  
    238282  // Default for matrices : Lines = Y , Columns = X
    239283  macoli_ = 0;
    240   marowi_ = 1;
     284  marowi_ = veceli_ = 1;
    241285  // Update OK
    242286  ndim_ = ndim;
    … …  
    278322  macoli_ = a.macoli_;
    279323  marowi_ = a.marowi_;
     324  veceli_ = a.veceli_;
    280325  // Update OK
    281326  ndim_ = a.ndim_;
    … …  
    284329
    285330
    286 void BaseArray::SubArray(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step)
    287 {
    288   if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("BaseArray::SubArray( ... ) NDim Error") );
     331void BaseArray::UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const
     332{
     333  if ( (ndim > ndim_) || (ndim < 1) )
     334    throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) NDim Error") );
    289335  int k;
    290336  for(k=0; k<ndim; k++)
    291337    if ( (siz[k]*step[k]+pos[k]) > size_[k] ) 
    292       throw(SzMismatchError("BaseArray::SubArray( ... ) Size/Pos Error") );
     338      throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) Size/Pos Error") );
    293339  uint_8 offset = offset_;
    294340  for(k=0; k<ndim_; k++) {
    … …  
    296342    step[k] *= step_[k];
    297343  }
    298   string exm = "BaseArray::SubArray() ";
     344  string exm = "BaseArray::UpdateSubArraySizes() ";
    299345  if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
    300346     throw( ParmError(exm) );
    301   //  (ra.mNDBlock).Share(mNDBlock);
    302347  return;
    303348}

  • trunk/SophyaLib/TArray/basarr.h

    r787 r804  
    2121class BaseArray : public AnyDataObj {
    2222public:
     23  //  To define Array / Matrix memory mapping
     24  enum { AutoMemoryMapping = -1, SameMemoryMapping = 0,
     25         CMemoryMapping = 1, FortranMemoryMapping = 2 };
     26
     27  //    Max Nb of printed elements
     28  static void    SetMaxPrint(uint_4 nprt=50);   
     29  //   Memory organisation (for matrices)
     30  static short   SetDefaultMemoryMapping(short mm=CMemoryMapping);
     31  static inline short GetDefaultMemoryMapping() { return default_memory_mapping; }
     32 
    2333  // Creation / destruction
    2434  BaseArray();
    … …  
    4454
    4555  // memory organisation - packing information
     56  inline short  GetMemoryMapping() const
     57                {return ( (marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping) ; }
     58
     59  inline uint_4 RowsKA() const {return marowi_; }    // Dimension correspondant aux lignes
     60  inline uint_4 ColsKA() const {return macoli_; }    // Dimension correspondant aux colonnes
     61  inline uint_4 VectKA() const {return veceli_; }    // Dimension corr. aux elts d'un vecteur
     62
    4663  inline bool   IsPacked() const { return(moystep_ == 1); }
    4764  inline bool   IsPackedX() const { return(step_[0] == 1); }
    … …  
    6683
    6784// Impression, I/O, ...
    68   static void   SetMaxPrint(uint_4 nprt=50);   
    69   void          Show(ostream& os, bool si=false) const;
    70   inline void   Show() const { Show(cout); }
     85  void           Show(ostream& os, bool si=false) const;
     86  inline void    Show() const { Show(cout); }
    7187  virtual string DataType() const = 0;
    7288
    … …  
    8399  bool UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg);
    84100  bool UpdateSizes(const BaseArray& a, string & exmsg);
    85   static uint_8 ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset) ;
    86   // Pour Extraction de sous-tableau
    87   virtual void SubArray(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step);
     101  static uint_8 ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset) ;
     102  //  Organisation memoire
     103  static short SelectMemoryMapping(short mm);
     104  void UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny);
     105  void UpdateMemoryMapping(BaseArray const & a, short mm);
     106
     107  // Pour Extraction de sous-tableau
     108  virtual void UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const;
    88109
    89110  uint_4 ndim_;                   // nb of dimensions
    … …  
    95116  uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
    96117  uint_4 marowi_, macoli_;     // For matrices, Row index and column index in dimensions
     118  uint_4 veceli_;              // For vectors,  dimension index = marowi_/macoli_ (Row/Col vectors)
    97119
    98120  DVList* mInfo;               // Infos (variables) attachees au tableau
    … …  
    100122  static char * ck_op_msg_[6];  // Operation messages for CheckDI() CheckBound()
    101123  static uint_4 max_nprt_;      // Nb maxi d'elements imprimes
     124  static short default_memory_mapping;  // Default memory mapping
    102125};
    103126

  • trunk/SophyaLib/TArray/fioarr.cc

    r787 r804  
    33
    44#include "pexceptions.h"
     5#include "fiondblock.h"
    56#include "fioarr.h"
     7#include "tmatrix.h"
     8#include "tvector.h"
    69
    710// --------------------------------------------------------
    … …  
    1316FIO_TArray<T>::FIO_TArray()
    1417{
    15   dobj=new TArray<T>;
    16   ownobj=true;
     18  dobj=NULL;
     19  ownobj=false;
    1720}
    1821
    … …  
    2023FIO_TArray<T>::FIO_TArray(string const & filename)
    2124{
    22   dobj=new TArray<T>;
    23   ownobj=true;
     25  dobj=NULL;
     26  ownobj=false;
    2427  Read(filename);
    2528}
    … …  
    2831FIO_TArray<T>::FIO_TArray(const TArray<T> & obj)
    2932{
    30   dobj = new TArray<T>(obj);
     33  const TVector<T> * tv = dynamic_cast<const TVector<T> *>(&obj);
     34  if (tv != NULL)  dobj = new TVector<T>(*tv, true);
     35  else {
     36  const TMatrix<T> * tm = dynamic_cast<const TMatrix<T> *>(&obj);
     37  if (tm != NULL)  dobj = new TMatrix<T>(*tm, true);
     38  else dobj = new TArray<T>(obj, true);
     39  }
    3140  ownobj=true;
    3241}
    … …  
    6372void FIO_TArray<T>::ReadSelf(PInPersist& is)
    6473{
    65 // On lit les 3 premiers uint_4
    66 //  0: Numero de version,   : NRows,  2 : NCol
     74// On lit les 5 premiers uint_4
     75//  0: Numero de version,  1 : Type (Array, matrix, Vector, ...)  2 != 0 , has Info
     76//          1:Type = 0 TArray , 12=(4+8) TMatrix , 48=(16+32) TVector
    6777  uint_4 itab[5];
    6878  is.Get(itab,5);
    69   if (dobj == NULL) dobj = new TArray<T>;
     79  if (dobj == NULL) {
     80    if (itab[1] == 12)   dobj = new TMatrix<T>;
     81    else if (itab[1] == 48)   dobj = new TVector<T>;
     82    else   dobj = new TArray<T>;
     83  }
    7084// On lit les tailles, etc ...
    7185  is.Get(dobj->ndim_);
    … …  
    7892  is.Get(dobj->marowi_);
    7993  is.Get(dobj->macoli_);
     94  is.Get(dobj->veceli_);
    8095// On lit le datablock
    8196  is >> dobj->DataBlock();
    … …  
    88103{
    89104  if (dobj == NULL)   return;
    90 //  On ecrit 4 uint_4 ....
     105//  On ecrit 5 uint_4 ....
    91106//  0: Numero de version,  1 : Type (Array, matrix, Vector, ...)  2 != 0 , has Info
     107//          1:Type = 0 TArray , 12=(4+8) TMatrix , 48=(16+32) TVector
     108  uint_4 typa = 0;
     109  TVector<T> * tv = dynamic_cast<TVector<T> *>(dobj);
     110  if (tv != NULL)  typa = 48;
     111  else {
     112  TMatrix<T> * tm = dynamic_cast<TMatrix<T> *>(dobj);
     113  if (tm != NULL)  typa = 12;
     114  else  typa = 0;
     115  }
     116
    92117  uint_4 itab[5];
    93118  itab[0] = 1;  // Numero de version a 1
    94   itab[1] = 0;
     119  itab[1] = typa;  // Real object type
    95120  itab[2] = (dobj->mInfo != NULL) ? 1 : 0;
    96121  itab[3] = itab[4] = 0;
    97   os.Put(itab,3);
     122  os.Put(itab,5);
    98123// On ecrit les tailles, etc ...
    99124  os.Put(dobj->ndim_);
    … …  
    106131  os.Put(dobj->marowi_);
    107132  os.Put(dobj->macoli_);
     133  os.Put(dobj->veceli_);
    108134// On ecrit le datablock
    109135  os << dobj->DataBlock();
    … …  
    117143#ifdef __CXX_PRAGMA_TEMPLATES__
    118144// Instances des delegues FileIO (PPersist)
    119 #pragma define_template FIO_TArray<uint_1>
     145// #pragma define_template FIO_TArray<uint_1>
    120146#pragma define_template FIO_TArray<uint_2>
    121 #pragma define_template FIO_TArray<int_2>
     147// #pragma define_template FIO_TArray<int_2>
    122148#pragma define_template FIO_TArray<int_4>
    123149#pragma define_template FIO_TArray<int_8>
    124 #pragma define_template FIO_TArray<uint_4>
    125 #pragma define_template FIO_TArray<uint_8>
     150// #pragma define_template FIO_TArray<uint_4>
     151// #pragma define_template FIO_TArray<uint_8>
    126152#pragma define_template FIO_TArray<r_8>
    127153#pragma define_template FIO_TArray<r_4>
    … …  
    132158#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    133159// Instances des delegues FileIO (PPersist)
    134 template class FIO_TArray<uint_1>;
     160// template class FIO_TArray<uint_1>;
    135161template class FIO_TArray<uint_2>;
    136 template class FIO_TArray<int_2>;
    137 template class FIO_TArray<int_4>;
     162// template class FIO_TArray<int_2>;
     163// template class FIO_TArray<int_4>;
    138164template class FIO_TArray<int_8>;
    139 template class FIO_TArray<uint_4>;
    140 template class FIO_TArray<uint_8>;
     165// template class FIO_TArray<uint_4>;
     166// template class FIO_TArray<uint_8>;
    141167template class FIO_TArray<r_8>;
    142168template class FIO_TArray<r_4>;

  • trunk/SophyaLib/TArray/fioarr.h

    r772 r804  
    99#include <iostream.h>
    1010#include "tarray.h"
     11#include "tmatrix.h"
     12#include "tvector.h"
    1113#include "ppersist.h"
    1214
    … …  
    3941{ FIO_TArray<T> fio(&obj);  fio.Read(is);  return(is); }
    4042
     43
    4144} // Fin du namespace
    4245

  • trunk/SophyaLib/TArray/matharr.cc

    r787 r804  
    1313TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_DoubleFunctionOfX f)
    1414{
     15  if (a.NbDimensions() < 1)
     16    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
    1517  T * pe;
    1618  uint_8 j,k;
    … …  
    3638TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_FloatFunctionOfX f)
    3739{
     40  if (a.NbDimensions() < 1)
     41    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
    3842  T * pe;
    3943  uint_8 j,k;
    … …  
    7579}
    7680
     81template <class T>
     82double MathArray<T>::MeanSigma(TArray<T> const & a, double & mean, double & sig)
     83{
     84  if (a.NbDimensions() < 1)
     85    throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
     86  const T * pe;
     87  uint_8 j,k;
     88  mean=0.;
     89  sig = 0.;
     90  double valok;
     91  if (a.AvgStep() > 0)   {  // regularly spaced elements
     92    uint_8 step = a.AvgStep();
     93    uint_8 maxx = a.Size()*step;
     94    pe = a.Data();
     95    for(k=0; k<maxx; k+=step )  {
     96      valok = (double) pe[k];
     97      mean += valok;  sig += valok*valok;
     98    }
     99  }
     100  else {    // Non regular data spacing ...
     101    uint_4 ka = a.MaxSizeKA();
     102    uint_8 step = a.Step(ka);
     103    uint_8 gpas = a.Size(ka)*step;
     104    for(j=0; j<a.Size(); j += a.Size(ka))  {
     105      pe = a.DataBlock().Begin()+a.Offset(j);
     106      for(k=0; k<gpas; k+=step) {
     107        valok = (double) pe[k];
     108        mean += valok;  sig += valok*valok;
     109      }
     110    }
     111  }
     112  double dsz = (double)(a.Size());
     113  mean /= dsz;
     114  sig = sig/dsz - mean*mean;
     115  if (sig >= 0.) sig = sqrt(sig);
     116  return(mean);
     117}
     118
    77119///////////////////////////////////////////////////////////////
    78120#ifdef __CXX_PRAGMA_TEMPLATES__

  • trunk/SophyaLib/TArray/matharr.h

    r787 r804  
    2525  virtual TArray<T>  ApplyFunction(TArray<T> const & a, Arr_DoubleFunctionOfX f);
    2626  virtual TArray<T>  ApplyFunction(TArray<T> const & a, Arr_FloatFunctionOfX f);
     27  // Computing Mean-Sigma
     28  virtual double     MeanSigma(TArray<T> const & a, double & mean, double & sig);
    2729};
    2830
    … …  
    4143  { MathArray<T> ma;   return( ma.ApplyFunction(a, tan) ); }
    4244
     45template <class T>
     46inline TArray<T> asin(const TArray<T>& a)
     47  { MathArray<T> ma;   return( ma.ApplyFunction(a, asin) ); }
     48
     49template <class T>
     50inline TArray<T> acos(const TArray<T>& a)
     51  { MathArray<T> ma;   return( ma.ApplyFunction(a, acos) ); }
     52
     53template <class T>
     54inline TArray<T> atan(const TArray<T>& a)
     55  { MathArray<T> ma;   return( ma.ApplyFunction(a, atan) ); }
     56
     57template <class T>
     58inline TArray<T> exp(const TArray<T>& a)
     59  { MathArray<T> ma;   return( ma.ApplyFunction(a, exp) ); }
     60
     61template <class T>
     62inline TArray<T> log(const TArray<T>& a)
     63  { MathArray<T> ma;   return( ma.ApplyFunction(a, log) ); }
     64
     65template <class T>
     66inline TArray<T> log10(const TArray<T>& a)
     67  { MathArray<T> ma;   return( ma.ApplyFunction(a, log10) ); }
     68
     69template <class T>
     70inline double MeanSigma(const TArray<T>& a, double & mean, double & sig)
     71  { MathArray<T> ma;   return( ma.MeanSigma(a, mean, sig) ); }
     72
     73template <class T>
     74inline double Mean(const TArray<T>& a)
     75  { MathArray<T> ma;  double mean, sig;  return( ma.MeanSigma(a, mean, sig) ); }
     76
     77
    4378} // Fin du namespace
    4479

  • trunk/SophyaLib/TArray/sopemtx.cc

    r772 r804  
    1414// -------------------------------------------------------------
    1515////////////////////////////////////////////////////////////////
     16/////////////////////////////////////////////////////////////////////////
     17// Classe de lignes/colonnes de matrices
     18enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
     19template <class T>
     20class TMatrixRC {
     21public:
     22  TMatrixRC();
     23  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
     24  virtual ~TMatrixRC() {}
     25
     26  // Acces aux rangees et colonnes de matrices
     27  static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
     28  static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
     29  static TMatrixRC<T> Diag(TMatrix<T> & m);
     30
     31  // ---- A virer   $CHECK$ Reza 03/2000
     32  //  int_4 Next();
     33  //  int_4 Prev();
     34  int_4 SetCol(int_4 c);
     35  int_4 SetRow(int_4 r);
     36  int_4 SetDiag();
     37
     38  static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
     39  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
     40
     41  TRCKind Kind() const { return kind; }
     42  uint_4 NElts() const;
     43  T& operator()(uint_4 i);
     44  T  operator()(uint_4 i) const;
     45
     46
     47   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
     48
     49// ---- A virer   $CHECK$ Reza 03/2000
     50//   TVector<T> GetVect() const;
     51//   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
     52//   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
     53
     54  TMatrixRC<T>& operator *= (T x);
     55  TMatrixRC<T>& operator /= (T x);
     56
     57//   TMatrixRC<T>& operator -= (T x);
     58//   TMatrixRC<T>& operator += (T x);
     59 
     60
     61  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
     62  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
     63
     64  uint_4 IMaxAbs(uint_4 first=0);
     65
     66  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
     67
     68  inline static double Abs_Value(uint_1 v) {return (double) v;}
     69  inline static double Abs_Value(uint_2 v) {return (double) v;}
     70  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
     71  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
     72  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
     73  inline static double Abs_Value(uint_4 v) {return (double) v;}
     74  inline static double Abs_Value(uint_8 v) {return (double) v;}
     75  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
     76  inline static double Abs_Value(r_8 v)    {return fabs(v);}
     77  inline static double Abs_Value(complex<float> v)
     78                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
     79  inline static double Abs_Value(complex<double> v)
     80                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
     81
     82protected:
     83  TMatrix<T>* matrix;
     84  T*          data;
     85  int_4       index;
     86  uint_4      step;
     87  TRCKind     kind;
     88};
     89
     90
     91
     92template <class T>
     93inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
     94  {
     95  if ( a.NElts() != b.NElts() )
     96    throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
     97  if ( a.Kind() != b.Kind() )
     98    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
     99  T sum = 0;
     100  for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
     101  return sum;
     102  }
     103
     104
     105template <class T>
     106inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
     107  { switch (rckind) { case TmatrixRow  : return m.Step(m.RowsKA());
     108                      case TmatrixCol  : return m.Step(m.ColsKA());
     109                      case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
     110    return 0; }
     111
     112template <class T>
     113inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
     114  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(&(m(index,0)));
     115                      case TmatrixCol  : return const_cast<T *>(&(m(0,index)));
     116                      case TmatrixDiag : return const_cast<T *>(m.Data()); }
     117    return NULL; }
     118
     119template <class T> inline uint_4 TMatrixRC<T>::NElts() const
     120  { if (!matrix) return 0;
     121    switch (kind) { case TmatrixRow  : return matrix->NCols();
     122                    case TmatrixCol  : return matrix->NRows();
     123                    case TmatrixDiag : return matrix->NCols(); }
     124    return 0; }
     125
     126template <class T>
     127inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
     128template <class T>
     129inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
     130
     131////////////////////////////////////////////////////////////////
     132// Typedef pour simplifier et compatibilite Peida
     133typedef TMatrixRC<r_8> MatrixRC;
     134
    16135
    17136template <class T> TMatrixRC<T>::TMatrixRC()
    … …  
    52171
    53172
    54 template <class T> int_4 TMatrixRC<T>::Next()
    55 {
    56 if (!matrix || kind==TmatrixDiag) return -1;
    57 index++;
    58 if(kind == TmatrixRow) {
    59   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
    60   data += matrix->NCols();
    61 } else {
    62   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
    63   data++;
    64 }
    65 return index;
    66 }
    67 
    68 template <class T> int_4 TMatrixRC<T>::Prev()
    69 {
    70 if (!matrix || kind == TmatrixDiag) return -1;
    71 index--;
    72 if(index < 0) {index = 0; return -1;}
    73 if(kind == TmatrixRow) data -= matrix->NCols();
    74 else data--;
    75 return index;
    76 }
    77 
     173// ---- A virer   $CHECK$ Reza 03/2000
     174// template <class T> int_4 TMatrixRC<T>::Next()
     175// {
     176// if (!matrix || kind==TmatrixDiag) return -1;
     177// index++;
     178// if(kind == TmatrixRow) {
     179//   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
     180//   data += matrix->NCols();
     181// } else {
     182//   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
     183//   data++;
     184// }
     185// return index;
     186// }
     187
     188// template <class T> int_4 TMatrixRC<T>::Prev()
     189// {
     190// if (!matrix || kind == TmatrixDiag) return -1;
     191// index--;
     192// if(index < 0) {index = 0; return -1;}
     193// if(kind == TmatrixRow) data -= matrix->NCols();
     194// else data--;
     195// return index;
     196// }
     197
     198 
    78199template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
    79200{
    … …  
    121242}
    122243
    123 template <class T> TVector<T> TMatrixRC<T>::GetVect() const
    124 {
    125 TVector<T> v(NElts());
    126 for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
    127 return v;
    128 }
    129 
    130 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
    131 {
    132 if ( NElts() != rc.NElts() )
    133   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
    134 if ( kind != rc.kind )
    135   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
    136 for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
    137 return *this;
    138 }
    139 
    140 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
    141 {
    142 if( NElts() != rc.NElts() )
    143   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
    144 if( kind != rc.kind )
    145   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
    146 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
    147 return *this;
    148 }
     244// ---- A virer   $CHECK$ Reza 03/2000
     245// template <class T> TVector<T> TMatrixRC<T>::GetVect() const
     246// {
     247// TVector<T> v(NElts());
     248// for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
     249// return v;
     250// }
     251
     252// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
     253// {
     254// if ( NElts() != rc.NElts() )
     255//   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
     256// if ( kind != rc.kind )
     257//   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
     258// for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
     259// return *this;
     260// }
     261
     262// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
     263// {
     264// if( NElts() != rc.NElts() )
     265//   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
     266// if( kind != rc.kind )
     267//   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
     268// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
     269// return *this;
     270// }
    149271
    150272
    … …  
    162284
    163285
    164 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
    165 {
    166 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
    167 return *this;
    168 }
    169 
    170 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
    171 {
    172 for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
    173 return *this;
    174 }
     286// ---- A virer   $CHECK$ Reza 03/2000
     287// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
     288// {
     289// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
     290// return *this;
     291// }
     292
     293// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
     294// {
     295// for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
     296// return *this;
     297// }
    175298
    176299template <class T>
    … …  
    331454
    332455
     456LinFitter::LinFitter()
     457{
     458}
     459
     460LinFitter::~LinFitter()
     461{
     462}
     463
     464double LinFitter::LinFit(const Vector& x, const Vector& y, int nf,
     465                         double (*f)(int, double), Vector& c)
     466{
     467  int n = x.NElts();
     468  if (n != y.NElts()) THROW(sizeMismatchErr);
     469 
     470  Matrix fx(nf, n);
     471  for (int i=0; i<nf; i++)
     472    for (int j=0; j<n; j++)
     473      fx(i,j) = f(i,x(j));
     474
     475  return LinFit(fx,y,c);
     476}
     477
     478
     479
     480double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
     481{
     482  int n = y.NElts();
     483  if ( n != fx.NCol()) THROW(sizeMismatchErr);
     484
     485  int nf = fx.NRows();
     486
     487  Matrix a(nf,nf);
     488
     489  for (int j=0; j<nf; j++)
     490    for (int k=j; k<nf; k++)
     491      a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j)
     492
     493        * TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
     494   
     495  c = fx * y;
     496
     497  if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
     498
     499  double xi2 = 0;
     500
     501  for (int k=0; k<n; k++) {
     502    double x = 0;
     503    for (int i=0; i<nf; i++)
     504      x += c(i) * fx(i,k);
     505    x -= y(k);
     506    xi2 += x*x;
     507  }
     508  return xi2;
     509}
     510
     511
     512
     513double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
     514                         double (*f)(int, double), Vector& c, Vector& errC)
     515{
     516  int n = x.NElts();
     517  if (n != y.NElts()) THROW(sizeMismatchErr);
     518 
     519  Matrix fx(nf, n);
     520  for (int i=0; i<nf; i++)
     521    for (int j=0; j<n; j++)
     522      fx(i,j) = f(i,x(j));
     523
     524  return LinFit(fx,y,errY2,c,errC);
     525}
     526
     527
     528double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
     529           Vector& c, Vector& errC)
     530{
     531  int n = y.NElts();
     532  if ( n != errY2.NElts()) THROW(sizeMismatchErr);
     533  if ( n != fx.NCol()) THROW(sizeMismatchErr);
     534
     535  int nf = fx.NRows();
     536
     537  Matrix a(nf,nf);
     538
     539  c.Realloc(nf);
     540  errC.Realloc(nf);
     541
     542  for (int j=0; j<nf; j++)
     543    for (int k=j; k<nf; k++) {
     544      double x=0;
     545      for (int l=0; l<n; l++)
     546        x += fx(j,l) * fx(k,l) / errY2(l);        // Matrice a inverser
     547      a(j,k) = a(k,j) = x;
     548    }
     549 
     550  Matrix d(nf,nf+1);
     551  for (int k=0; k<nf; k++) {
     552    double x=0;
     553    for (int l=0; l<n; l++)
     554      x += fx(k,l) * y(l) / errY2(l);             // Second membre 1ere colonne
     555    d(k,0) = x;
     556    for (int m=1; m<=nf; m++)
     557      d(k,m) = (k==m) ? 1 : 0;                // Reste second membre = Id.
     558  }
     559
     560  if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
     561
     562  for (int l=0; l<nf; l++) {
     563    c(l) = d(l,0);                            // Parametres = 1ere colonne
     564    errC(l) = d(l,l+1);                        // Erreurs = diag inverse.
     565  }
     566
     567  double xi2 = 0;
     568
     569  for (int jj=0; jj<n; jj++) {
     570    double x = 0;
     571    for (int ii=0; ii<nf; ii++)
     572      x += c(ii) * fx(ii,jj);
     573    x -= y(jj);
     574    xi2 += x*x/errY2(jj);
     575  }
     576  return xi2;
     577}
     578
     579
     580
    333581
    334582///////////////////////////////////////////////////////////////
    335583#ifdef __CXX_PRAGMA_TEMPLATES__
    336584// Instances gestion lignes/colonnes
    337 #pragma define_template TMatrixRC<uint_1>
    338 #pragma define_template TMatrixRC<uint_2>
    339 #pragma define_template TMatrixRC<int_2>
    340585#pragma define_template TMatrixRC<int_4>
    341 #pragma define_template TMatrixRC<int_8>
    342 #pragma define_template TMatrixRC<uint_4>
    343 #pragma define_template TMatrixRC<uint_8>
    344586#pragma define_template TMatrixRC<r_4>
    345587#pragma define_template TMatrixRC<r_8>
    346 #pragma define_template TMatrixRC< complex<r_4> >
    347 #pragma define_template TMatrixRC< complex<r_8> >
    348 // #pragma define_template SimpleMatrixOperation<r_4>
     588#pragma define_template SimpleMatrixOperation<int_4>
     589#pragma define_template SimpleMatrixOperation<r_4>
    349590#pragma define_template SimpleMatrixOperation<r_8>
    350591#endif
    … …  
    352593#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    353594// Instances gestion lignes/colonnes
    354 template class TMatrixRC<uint_1>;
    355 template class TMatrixRC<uint_2>;
    356 template class TMatrixRC<int_2>;
    357595template class TMatrixRC<int_4>;
    358 template class TMatrixRC<int_8>;
    359 template class TMatrixRC<uint_4>;
    360 template class TMatrixRC<uint_8>;
    361596template class TMatrixRC<r_4>;
    362597template class TMatrixRC<r_8>;
    363 template class TMatrixRC< complex<r_4> >;
    364 template class TMatrixRC< complex<r_8> >;
     598template class SimpleMatrixOperation<int_4>;
    365599template class SimpleMatrixOperation<r_4>;
    366600template class SimpleMatrixOperation<r_8>;

  • trunk/SophyaLib/TArray/sopemtx.h

    r772 r804  
    99namespace SOPHYA {
    1010
    11 /////////////////////////////////////////////////////////////////////////
    12 // Classe de lignes/colonnes de matrices
    13 enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
    14 template <class T>
    15 class TMatrixRC {
    16 public:
    17   TMatrixRC();
    18   TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
    19   virtual ~TMatrixRC() {}
    2011
    21   // Acces aux rangees et colonnes de matrices
    22   static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
    23   static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
    24   static TMatrixRC<T> Diag(TMatrix<T> & m);
    25 
    26   int_4 Next();
    27   int_4 Prev();
    28   int_4 SetCol(int_4 c);
    29   int_4 SetRow(int_4 r);
    30   int_4 SetDiag();
    31 
    32   static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
    33   static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
    34 
    35   TRCKind Kind() const { return kind; }
    36   uint_4 NElts() const;
    37   T& operator()(uint_4 i);
    38   T  operator()(uint_4 i) const;
    39 
    40   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
    41   TVector<T> GetVect() const;
    42 
    43   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
    44   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
    45 
    46   TMatrixRC<T>& operator *= (T x);
    47   TMatrixRC<T>& operator /= (T x);
    48   TMatrixRC<T>& operator -= (T x);
    49   TMatrixRC<T>& operator += (T x);
    50 
    51   TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
    52   TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
    53 
    54   uint_4 IMaxAbs(uint_4 first=0);
    55 
    56   static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
    57 
    58   inline static double Abs_Value(uint_1 v) {return (double) v;}
    59   inline static double Abs_Value(uint_2 v) {return (double) v;}
    60   inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
    61   inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
    62   inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
    63   inline static double Abs_Value(uint_4 v) {return (double) v;}
    64   inline static double Abs_Value(uint_8 v) {return (double) v;}
    65   inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
    66   inline static double Abs_Value(r_8 v)    {return fabs(v);}
    67   inline static double Abs_Value(complex<float> v)
    68                 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    69   inline static double Abs_Value(complex<double> v)
    70                 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    71 
    72 protected:
    73   TMatrix<T>* matrix;
    74   T*          data;
    75   int_4       index;
    76   uint_4      step;
    77   TRCKind     kind;
    78 };
    79 
    80 
    81 template <class T>
    82 inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
    83   {
    84   if ( a.NElts() != b.NElts() )
    85     throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
    86   if ( a.Kind() != b.Kind() )
    87     throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
    88   T sum = 0;
    89   for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
    90   return sum;
    91   }
    92 
    93 template <class T>
    94 inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
    95   { switch (rckind) { case TmatrixRow  : return 1;
    96                       case TmatrixCol  : return m.NCols();
    97                       case TmatrixDiag : return m.NCols()+1; }
    98     return 0; }
    99 
    100 template <class T>
    101 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
    102   { switch (rckind) { case TmatrixRow  : return const_cast<T *>(m.Data()) + index * m.NCols();
    103                       case TmatrixCol  : return const_cast<T *>(m.Data()) + index;
    104                       case TmatrixDiag : return const_cast<T *>(m.Data()); }
    105     return NULL; }
    106 
    107 template <class T> inline uint_4 TMatrixRC<T>::NElts() const
    108   { if (!matrix) return 0;
    109     switch (kind) { case TmatrixRow  : return matrix->NCols();
    110                     case TmatrixCol  : return matrix->NRows();
    111                     case TmatrixDiag : return matrix->NCols(); }
    112     return 0; }
    113 
    114 template <class T>
    115 inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
    116 template <class T>
    117 inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
    118 
    119 ////////////////////////////////////////////////////////////////
    120 // Typedef pour simplifier et compatibilite Peida
    121 typedef TMatrixRC<r_8> MatrixRC;
    12212
    12313////////////////////////////////////////////////////////////////
    … …  
    14939}
    15040
     41//--------------------------------------
     42//        Linear fitting
     43//--------------------------------------
     44
     45class LinFitter {
     46public :
     47                LinFitter();
     48  virtual       ~LinFitter();
     49
     50  double        LinFit(const Vector& x, const Vector& y, int nf,
     51                       double (*f)(int, double), Vector& c);
     52// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1;
     53
     54  double        LinFit(const Matrix& fx, const Vector& y, Vector& c);
     55// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
     56// la matrice fx contient les valeurs des f:
     57// fx(i,j) = f(i, x(j)).
     58                     
     59  double        LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
     60                       double (*f)(int, double), Vector& c, Vector& errC);
     61// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
     62// errY2 contient les carres des erreurs sur les Y.
     63// au retour, errC contient les erreurs sur les coefs.
     64
     65  double        LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
     66                       Vector& c, Vector& errC);
     67// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
     68// la matrice fx contient les valeurs des f:
     69// fx(i,j) = f(i, x(j)).
     70// errY2 contient les carres des erreurs sur les Y.
     71// au retour, errC contient les erreurs sur les coefs.
     72};
     73
    15174
    15275} // Fin du namespace

  • trunk/SophyaLib/TArray/tarray.cc

    r787 r804  
    1818template <class T>
    1919TArray<T>::TArray()
    20   : mNDBlock() , BaseArray()
     20  : BaseArray() , mNDBlock()
    2121// Default constructor
    2222{
    … …  
    2424
    2525template <class T>
    26 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
    27   : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , BaseArray()
     26TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, uint_4 step)
     27  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1))
    2828{
    2929  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
    … …  
    3232
    3333template <class T>
    34 TArray<T>::TArray(uint_4 nx, uint_4 ny=1, uint_4 nz=1)
    35   : mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)) , BaseArray()
     34TArray<T>::TArray(uint_4 nx, uint_4 ny, uint_4 nz, uint_4 nt, uint_4 nu)
     35  : BaseArray() , mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)*((nt>0)?nt:1)*((nu>0)?nu:1))
    3636{
    3737  uint_4 size[BASEARRAY_MAXNDIMS];
    3838  size[0] = nx;   size[1] = ny;   size[2] = nz;
     39  size[3] = nt;   size[4] = nu;
    3940  int ndim = 1;
    40   if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
     41  if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) && (size[4] > 0) ) ndim = 5;
     42  else if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) ) ndim = 4;
     43  else if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
    4144  else if (size[1] > 0)  ndim = 2;
    4245  else ndim = 1;
    … …  
    4649
    4750template <class T>
    48 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
    49   : mNDBlock(db, share) , BaseArray()
     51TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
     52  : BaseArray() , mNDBlock(db, share)
    5053{
    5154  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
    … …  
    5558
    5659template <class T>
    57 TArray<T>::TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
    58   : mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) , BaseArray()
     60TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
     61  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br)
    5962{
    6063  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
    … …  
    6467template <class T>
    6568TArray<T>::TArray(const TArray<T>& a)
    66   : mNDBlock(a.mNDBlock) , BaseArray()
     69  : BaseArray() , mNDBlock(a.mNDBlock)
    6770{
    6871  string exmsg = "TArray<T>::TArray(const TArray<T>&)";
    … …  
    7376template <class T>
    7477TArray<T>::TArray(const TArray<T>& a, bool share)
    75   : mNDBlock(a.mNDBlock, share) , BaseArray()
     78  : BaseArray() , mNDBlock(a.mNDBlock, share)
    7679{
    7780  string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
    … …  
    8790
    8891template <class T>
    89 TArray<T>& TArray<T>::operator = (const TArray<T>& a)
     92TArray<T>& TArray<T>::Set(const TArray<T>& a)
    9093{
    9194  if (this != &a) {
    … …  
    171174
    172175
     176template <class T>
     177TArray<T> TArray<T>::PackElements(bool force) const
     178{
     179  if (NbDimensions() < 1)
     180    throw RangeCheckError("TArray<T>::PackElements() - Not Allocated Array ! ");
     181  if ( !force && (AvgStep() == 1) ) {
     182    TArray<T> ra(*this, true);
     183    ra.SetTemp(true);
     184    return(ra);
     185  }
     186  else {
     187    TArray<T> ra(ndim_, size_, 1);
     188    ra.CopyElt(*this);
     189    ra.SetTemp(true);
     190    return(ra);
     191  }
     192}
     193
    173194// SubArrays
    174 template <class T>
    175 TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
    176 {
     195// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     196template <class T>
     197TArray<T> TArray<T>::SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const
     198{
     199  if (NbDimensions() < 1)
     200    throw RangeCheckError("TArray<T>::operator () (Range, ...) - Not Allocated Array ! ");
    177201  uint_4 ndim = 0;
    178202  uint_4 size[BASEARRAY_MAXNDIMS];
    … …  
    199223  ndim = ndim_;
    200224  TArray<T> ra;
    201   SubArray(ra, ndim, size, pos, step);
     225  UpdateSubArraySizes(ra, ndim, size, pos, step);
    202226  ra.DataBlock().Share(this->DataBlock());
    203227  ra.SetTemp(true);
    … …  
    210234//  Possibilite de // , vectorisation
    211235template <class T>
    212 TArray<T>& TArray<T>::operator = (Sequence seq)
    213 {
     236TArray<T>& TArray<T>::Set(Sequence seq)
     237{
     238  if (NbDimensions() < 1)
     239    throw RangeCheckError("TArray<T>::Set(Sequence ) - Not Allocated Array ! ");
    214240  T * pe;
    215241  uint_8 j,k;
    … …  
    234260
    235261template <class T>
    236 TArray<T>& TArray<T>::operator = (T x)
    237 {
     262TArray<T>& TArray<T>::Set(T x)
     263{
     264  if (NbDimensions() < 1)
     265    throw RangeCheckError("TArray<T>::Set(T )  - Not Allocated Array ! ");
    238266  T * pe;
    239267  uint_8 j,k;
    … …  
    257285
    258286template <class T>
    259 TArray<T>& TArray<T>::operator += (T x)
    260 {
     287TArray<T>& TArray<T>::Add(T x)
     288{
     289  if (NbDimensions() < 1)
     290    throw RangeCheckError("TArray<T>::Add(T )  - Not Allocated Array ! ");
    261291  T * pe;
    262292  uint_8 j,k;
    … …  
    280310
    281311template <class T>
    282 TArray<T>& TArray<T>::operator -= (T x)
    283 {
     312TArray<T>& TArray<T>::Sub(T x)
     313{
     314  if (NbDimensions() < 1)
     315    throw RangeCheckError("TArray<T>::Sub(T )  - Not Allocated Array ! ");
    284316  T * pe;
    285317  uint_8 j,k;
    … …  
    303335
    304336template <class T>
    305 TArray<T>& TArray<T>::operator *= (T x)
    306 {
     337TArray<T>& TArray<T>::Mul(T x)
     338{
     339  if (NbDimensions() < 1)
     340    throw RangeCheckError("TArray<T>::Mul(T )  - Not Allocated Array ! ");
    307341  T * pe;
    308342  uint_8 j,k;
    … …  
    326360
    327361template <class T>
    328 TArray<T>& TArray<T>::operator /= (T x)
    329 {
     362TArray<T>& TArray<T>::Div(T x)
     363{
     364  if (NbDimensions() < 1)
     365    throw RangeCheckError("TArray<T>::Div(T )  - Not Allocated Array ! ");
    330366  T * pe;
    331367  uint_8 j,k;
    … …  
    349385
    350386
     387template <class T>
     388TArray<T>& TArray<T>::SubInv(T x)
     389{
     390  if (NbDimensions() < 1)
     391    throw RangeCheckError("TArray<T>::SubInv(T )  - Not Allocated Array ! ");
     392  T * pe;
     393  uint_8 j,k;
     394  if (AvgStep() > 0)   {  // regularly spaced elements
     395    uint_8 step = AvgStep();
     396    uint_8 maxx = totsize_*step;
     397    pe = Data();
     398    for(k=0; k<maxx; k+=step )  pe[k] = x-pe[k];
     399  }
     400  else {    // Non regular data spacing ...
     401    uint_4 ka = MaxSizeKA();
     402    uint_8 step = Step(ka);
     403    uint_8 gpas = Size(ka)*step;
     404    for(j=0; j<totsize_; j += Size(ka))  {
     405      pe = mNDBlock.Begin()+Offset(j);
     406      for(k=0; k<gpas; k+=step)  pe[k] = x-pe[k];
     407    }
     408  }
     409  return(*this);
     410}
     411
     412template <class T>
     413TArray<T>& TArray<T>::DivInv(T x)
     414{
     415  if (NbDimensions() < 1)
     416    throw RangeCheckError("TArray<T>::DivInv(T )  - Not Allocated Array ! ");
     417  T * pe;
     418  uint_8 j,k;
     419  if (AvgStep() > 0)   {  // regularly spaced elements
     420    uint_8 step = AvgStep();
     421    uint_8 maxx = totsize_*step;
     422    pe = Data();
     423    for(k=0; k<maxx; k+=step )  pe[k] = x/pe[k];
     424  }
     425  else {    // Non regular data spacing ...
     426    uint_4 ka = MaxSizeKA();
     427    uint_8 step = Step(ka);
     428    uint_8 gpas = Size(ka)*step;
     429    for(j=0; j<totsize_; j += Size(ka))  {
     430      pe = mNDBlock.Begin()+Offset(j);
     431      for(k=0; k<gpas; k+=step)  pe[k] = x/pe[k];
     432    }
     433  }
     434  return(*this);
     435}
     436
     437
    351438//  >>>> Operations avec 2nd membre de type tableau
    352439template <class T>
    353 TArray<T>& TArray<T>::operator += (const TArray<T>& a)
    354 {
    355   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator += (const TArray<T>&)")) ;
     440TArray<T>& TArray<T>::AddElt(const TArray<T>& a)
     441{
     442  if (NbDimensions() < 1)
     443    throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
     444  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&)")) ;
    356445
    357446  T * pe;
    … …  
    381470
    382471template <class T>
    383 TArray<T>& TArray<T>::operator -= (const TArray<T>& a)
    384 {
    385   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
     472TArray<T>& TArray<T>::SubElt(const TArray<T>& a)
     473{
     474  if (NbDimensions() < 1)
     475    throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
     476  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&)")) ;
    386477
    387478  T * pe;
    … …  
    412503
    413504template <class T>
    414 TArray<T>& TArray<T>::MultElt(const TArray<T>& a)
    415 {
    416   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
     505TArray<T>& TArray<T>::MulElt(const TArray<T>& a)
     506{
     507  if (NbDimensions() < 1)
     508    throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
     509  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&)")) ;
    417510
    418511  T * pe;
    … …  
    441534}
    442535
     536
    443537template <class T>
    444538TArray<T>& TArray<T>::DivElt(const TArray<T>& a)
    445539{
     540  if (NbDimensions() < 1)
     541    throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
    446542  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
    447543
    … …  
    471567}
    472568
     569template <class T>
     570TArray<T>& TArray<T>::CopyElt(const TArray<T>& a)
     571{
     572  if (NbDimensions() < 1)
     573    throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
     574  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
     575
     576  T * pe;
     577  const T * pea;
     578  uint_8 j,k,ka;
     579  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     580    uint_8 step = AvgStep();
     581    uint_8 stepa = a.AvgStep();
     582    uint_8 maxx = totsize_*step;
     583    pe = Data();
     584    pea = a.Data();
     585    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka] ;
     586  }
     587  else {    // Non regular data spacing ...
     588    uint_4 ax = MaxSizeKA();
     589    uint_8 step = Step(ax);
     590    uint_8 stepa = a.Step(ax);
     591    uint_8 gpas = Size(ax)*step;
     592    for(j=0; j<totsize_; j += Size(ax))  {
     593      pe = mNDBlock.Begin()+Offset(j);
     594      pea = a.DataBlock().Begin()+a.Offset(j);
     595      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka];
     596    }
     597  }
     598  return(*this);
     599}
     600
     601
     602// Somme et produit des elements
     603template <class T>
     604T TArray<T>::Sum() const
     605{
     606  if (NbDimensions() < 1)
     607    throw RangeCheckError("TArray<T>::Sum()  - Not Allocated Array ! ");
     608  T ret=0;
     609  const T * pe;
     610  uint_8 j,k;
     611  if (AvgStep() > 0)   {  // regularly spaced elements
     612    uint_8 step = AvgStep();
     613    uint_8 maxx = totsize_*step;
     614    pe = Data();
     615    for(k=0; k<maxx; k+=step )  ret += pe[k];
     616  }
     617  else {    // Non regular data spacing ...
     618    uint_4 ka = MaxSizeKA();
     619    uint_8 step = Step(ka);
     620    uint_8 gpas = Size(ka)*step;
     621    for(j=0; j<totsize_; j += Size(ka))  {
     622      pe = mNDBlock.Begin()+Offset(j);
     623      for(k=0; k<gpas; k+=step)  ret += pe[k] ;
     624    }
     625  }
     626  return ret;
     627}
     628
     629template <class T>
     630T TArray<T>::Product() const
     631{
     632  if (NbDimensions() < 1)
     633    throw RangeCheckError("TArray<T>::Product()  - Not Allocated Array ! ");
     634  T ret=0;
     635  const T * pe;
     636  uint_8 j,k;
     637  if (AvgStep() > 0)   {  // regularly spaced elements
     638    uint_8 step = AvgStep();
     639    uint_8 maxx = totsize_*step;
     640    pe = Data();
     641    for(k=0; k<maxx; k+=step )  ret *= pe[k];
     642  }
     643  else {    // Non regular data spacing ...
     644    uint_4 ka = MaxSizeKA();
     645    uint_8 step = Step(ka);
     646    uint_8 gpas = Size(ka)*step;
     647    for(j=0; j<totsize_; j += Size(ka))  {
     648      pe = mNDBlock.Begin()+Offset(j);
     649      for(k=0; k<gpas; k+=step)  ret *= pe[k] ;
     650    }
     651  }
     652  return ret;
     653}
     654
     655
    473656
    474657// ----------------------------------------------------
    … …  
    487670{
    488671  if (maxprt < 0)  maxprt = max_nprt_;
    489   uint_4 npr = 0;
     672  int_4 npr = 0;
    490673  Show(os, si);
    491   int k0,k1,k2,k3,k4;
     674  uint_4 k0,k1,k2,k3,k4;
    492675  for(k4=0; k4<size_[4]; k4++) {
    493676    if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
    … …  
    535718///////////////////////////////////////////////////////////////
    536719#ifdef __CXX_PRAGMA_TEMPLATES__
     720/*
    537721#pragma define_template TArray<uint_1>
     722#pragma define_template TArray<int_2>
     723#pragma define_template TArray<uint_4>
     724#pragma define_template TArray<uint_8>
     725*/
    538726#pragma define_template TArray<uint_2>
    539 #pragma define_template TArray<int_2>
    540727#pragma define_template TArray<int_4>
    541728#pragma define_template TArray<int_8>
    542 #pragma define_template TArray<uint_4>
    543 #pragma define_template TArray<uint_8>
    544729#pragma define_template TArray<r_4>
    545730#pragma define_template TArray<r_8>
    … …  
    549734
    550735#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    551 template class TArray<uint_1>;
     736/*
     737template class TArray<uint_1>;
     738template class TArray<int_2>; 
     739template class TArray<uint_4>;
     740template class TArray<uint_8>;
     741*/
    552742template class TArray<uint_2>;
    553 template class TArray<int_2>;
    554743template class TArray<int_4>;
    555744template class TArray<int_8>;
    556 template class TArray<uint_4>;
    557 template class TArray<uint_8>;
    558745template class TArray<r_4>;
    559746template class TArray<r_8>;

  • trunk/SophyaLib/TArray/tarray.h

    r787 r804  
    2828  // Creation / destruction
    2929  TArray();
    30   TArray(uint_4 ndim, uint_4 * siz, uint_4 step =1);
    31   TArray(uint_4 nx, uint_4 ny=0, uint_4 nz=0);
    32   TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share=false, uint_4 step=1, uint_8 offset=0);
    33   TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step=1, uint_8 offset=0, Bridge* br=NULL);
     30  TArray(uint_4 ndim, const uint_4 * siz, uint_4 step =1);
     31  TArray(uint_4 nx, uint_4 ny=0, uint_4 nz=0, uint_4 nt=0, uint_4 nu=0);
     32  TArray(uint_4 ndim, const uint_4 * siz, NDataBlock<T> & db, bool share=false, uint_4 step=1, uint_8 offset=0);
     33  TArray(uint_4 ndim, const uint_4 * siz, T* values, uint_4 step=1, uint_8 offset=0, Bridge* br=NULL);
    3434  TArray(const TArray<T>& a);
    3535  TArray(const TArray<T>& a, bool share);
    … …  
    3838
    3939  // A = B : partage les donnees si "a" est temporaire, clone sinon.
    40   virtual TArray<T>& operator = (const TArray<T>& a);
     40  inline  TArray<T>& operator = (const TArray<T>& a) { return Set(a); }
     41  virtual TArray<T>& Set(const TArray<T>& a);
    4142
    4243  // Gestion taille/Remplissage
    4344  virtual void Clone(const TArray<T>& a);
    44   virtual void ReSize(uint_4 ndim, uint_4 * siz, uint_4 step=1);
    45   virtual void Realloc(uint_4 ndim, uint_4 * siz, uint_4 step=1, bool force=false);
     45  void ReSize(uint_4 ndim, uint_4 * siz, uint_4 step=1);
     46  void Realloc(uint_4 ndim, uint_4 * siz, uint_4 step=1, bool force=false);
    4647
    4748  // Compacts size=1 array dimensions
    … …  
    4950  virtual TArray<T>& CompactTrailingDimensions();
    5051
    51   // SubArrays
    52   virtual TArray<T> operator () (Range rx, Range ry, Range rz=0, Range rt=0, Range ru=0);
    53 
     52  // Packing array elements in memory
     53  virtual TArray<T> PackElements(bool force=false) const ;
     54
     55  // SubArrays - $CHECK$ Reza 03/2000 je ne sais pas s'il faut declarer ca const ??
     56  TArray<T> SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const ;
     57  inline TArray<T> operator () (Range rx, Range ry, Range rz, Range rt=0, Range ru=0) const
     58                   { return  SubArray(rx, ry, rz, rt, ru); }
    5459
    5560  // ---- Access to data
    … …  
    8388  inline T toScalar();
    8489// Met les elements a une suite de valeurs
    85   virtual TArray<T>&  operator = (Sequence seq);
     90  virtual TArray<T>&  Set(Sequence seq);
     91  inline  TArray<T>&  operator = (Sequence seq)    { return Set(seq); }
    8692// A = x (tous les elements a x)
    87   virtual TArray<T>&  operator = (T x);
     93  virtual TArray<T>&  Set(T x);
     94  inline  TArray<T>&  operator = (T x)             { return Set(x); }
    8895// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
    89   virtual TArray<T>&  operator += (T x);
    90   virtual TArray<T>&  operator -= (T x);
    91   virtual TArray<T>&  operator *= (T x);
    92   virtual TArray<T>&  operator /= (T x);
     96  virtual TArray<T>&  Add(T x);
     97  inline  TArray<T>&  operator += (T x)            { return Add(x); }
     98  virtual TArray<T>&  Sub(T x);
     99  inline  TArray<T>&  operator -= (T x)            { return Sub(x); }
     100  virtual TArray<T>&  Mul(T x);
     101  inline  TArray<T>&  operator *= (T x)            { return Mul(x); }
     102  virtual TArray<T>&  Div(T x);
     103  inline  TArray<T>&  operator /= (T x)            { return Div(x); }
     104  virtual TArray<T>&  SubInv(T x);  //   A ---> x-A
     105  virtual TArray<T>&  DivInv(T x);  //   A ---> x/A
     106
    93107// A += -=  (ajoute, soustrait element par element les deux tableaux )
    94   virtual TArray<T>&  operator += (const TArray<T>& a);
    95   virtual TArray<T>&  operator -= (const TArray<T>& a);
     108  virtual TArray<T>&  AddElt(const TArray<T>& a);
     109  inline  TArray<T>&  operator += (const TArray<T>& a)  { return AddElt(a); }
     110  virtual TArray<T>&  SubElt(const TArray<T>& a);
     111  inline  TArray<T>&  operator -= (const TArray<T>& a)  { return SubElt(a); }
    96112// Multiplication, division element par element les deux tableaux
    97   virtual TArray<T>&  MultElt(const TArray<T>& a);
     113  virtual TArray<T>&  MulElt(const TArray<T>& a);
    98114  virtual TArray<T>&  DivElt(const TArray<T>& a);
     115// Recopie des valeurs, element par element
     116  virtual TArray<T>&  CopyElt(const TArray<T>& a);
     117
     118// Somme et produit des elements
     119  virtual T Sum() const ;
     120  virtual T Product() const ;
    99121
    100122// Impression, I/O, ...
    … …  
    124146
    125147template <class T> inline TArray<T> operator + (const TArray<T>& a, T b)
    126     {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
     148    {TArray<T> result(a); result.SetTemp(true); result.Add(b); return result;}
    127149
    128150template <class T> inline TArray<T> operator + (T b,const TArray<T>& a)
    129     {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
     151    {TArray<T> result(a); result.SetTemp(true); result.Add(b); return result;}
    130152
    131153template <class T> inline TArray<T> operator - (const TArray<T>& a, T b)
    132     {TArray<T> result(a); result.SetTemp(true); result -= b; return result;}
     154    {TArray<T> result(a); result.SetTemp(true); result.Sub(b); return result;}
    133155
    134156template <class T> inline TArray<T> operator - (T b,const TArray<T>& a)
    135     {TArray<T> result(a); result.SetTemp(true);
    136      result.DataBlock() = b-result.DataBlock(); return result;}
     157    {TArray<T> result(a); result.SetTemp(true); result.SubInv(b); return result;}
    137158
    138159template <class T> inline TArray<T> operator * (const TArray<T>& a, T b)
    139     {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}
     160    {TArray<T> result(a); result.SetTemp(true); result.Mul(b); return result;}
    140161
    141162template <class T> inline TArray<T> operator * (T b,const TArray<T>& a)
    142     {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}
     163    {TArray<T> result(a); result.SetTemp(true); result.Mul(b); return result;}
    143164
    144165template <class T> inline TArray<T> operator / (const TArray<T>& a, T b)
    145     {TArray<T> result(a); result.SetTemp(true); result /= b; return result;}
     166    {TArray<T> result(a); result.SetTemp(true); result.DivInv(b); return result;}
    146167
    147168////////////////////////////////////////////////////////////////
    … …  
    150171template <class T>
    151172inline TArray<T> operator + (const TArray<T>& a,const TArray<T>& b)
    152     {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
     173    {TArray<T> result(a); result.SetTemp(true); result.AddElt(b); return result;}
    153174
    154175template <class T>
    155176inline TArray<T> operator - (const TArray<T>& a,const TArray<T>& b)
    156     {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
     177    {TArray<T> result(a); result.SetTemp(true); result.SubElt(b); return result;}
    157178
    158179
    … …  
    180201{
    181202  CheckBound(ix, iy, iz, it, iu, 4);
    182   Elem(ix, iy, iz, it, iu);
     203  return(Elem(ix, iy, iz, it, iu));
    183204}
    184205
    … …  
    187208{
    188209  CheckBound(ix, iy, iz, it, iu, 4);
    189   Elem(ix, iy, iz, it, iu);
     210  return(Elem(ix, iy, iz, it, iu));
    190211}
    191212
    … …  
    259280}
    260281
     282// Typedef pour simplifier
     283typedef TArray<r_8> Array;
     284
    261285} // Fin du namespace
    262286

  • trunk/SophyaLib/TArray/tarrinit.cc

    r772 r804  
    33#include "tarrinit.h"
    44
    5 #include "tmatrix.h"
    6 #include "tvector.h"
    75#include "fioarr.h"
    86
    … …  
    1917 
    2018
    21   PPRegister(FIO_TArray<uint_1>);
     19  //  PPRegister(FIO_TArray<uint_1>);
     20  //  DObjRegister(FIO_TArray<uint_1>, TArray<uint_1>);
    2221  PPRegister(FIO_TArray<uint_2>);
    23   PPRegister(FIO_TArray<int_2>);
     22  DObjRegister(FIO_TArray<uint_2>, TArray<uint_2>);
     23  //  PPRegister(FIO_TArray<int_2>);
     24  //  DObjRegister(FIO_TArray<int_2>, TArray<int_2>);
    2425  PPRegister(FIO_TArray<int_4>);
     26  DObjRegister(FIO_TArray<int_4>, TArray<int_4>);
    2527  PPRegister(FIO_TArray<int_8>);
    26   PPRegister(FIO_TArray<uint_8>);
     28  DObjRegister(FIO_TArray<int_8>, TArray<int_8>);
     29  //  PPRegister(FIO_TArray<uint_4>);
     30  // DObjRegister(FIO_TArray<uint_4>, TArray<uint_4>);
     31  //  PPRegister(FIO_TArray<uint_8>);
     32  //  DObjRegister(FIO_TArray<uint_8>, TArray<uint_8>);
    2733  PPRegister(FIO_TArray<r_4>);
     34  DObjRegister(FIO_TArray<r_4>, TArray<r_4>);
    2835  PPRegister(FIO_TArray<r_8>);
     36  DObjRegister(FIO_TArray<r_8>, TArray<r_8>);
    2937  PPRegister(FIO_TArray< complex<r_4> >);
     38  DObjRegister(FIO_TArray< complex<r_4> >, TArray< complex<r_4> >);
    3039  PPRegister(FIO_TArray< complex<r_8> >);
    31 
    32   PPRegister(FIO_TMatrix<uint_1>);
    33   PPRegister(FIO_TMatrix<uint_2>);
    34   PPRegister(FIO_TMatrix<int_2>);
    35   PPRegister(FIO_TMatrix<int_4>);
    36   PPRegister(FIO_TMatrix<int_8>);
    37   PPRegister(FIO_TMatrix<uint_8>);
    38   PPRegister(FIO_TMatrix<r_4>);
    39   PPRegister(FIO_TMatrix<r_8>);
    40   PPRegister(FIO_TMatrix< complex<r_4> >);
    41   PPRegister(FIO_TMatrix< complex<r_8> >);
    42 
    43   PPRegister(FIO_TVector<uint_1>);
    44   PPRegister(FIO_TVector<uint_2>);
    45   PPRegister(FIO_TVector<int_2>);
    46   PPRegister(FIO_TVector<int_4>);
    47   PPRegister(FIO_TVector<int_8>);
    48   PPRegister(FIO_TVector<uint_8>);
    49   PPRegister(FIO_TVector<r_4>);
    50   PPRegister(FIO_TVector<r_8>);
    51   PPRegister(FIO_TVector< complex<r_4> >);
    52   PPRegister(FIO_TVector< complex<r_8> >);
    53 
    54 
     40  DObjRegister(FIO_TArray< complex<r_8> >, TArray< complex<r_8> >);
    5541
    5642}

  • trunk/SophyaLib/TArray/tmatrix.cc

    r772 r804  
    1 // $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
     1// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    44#include <stdio.h>
    55#include <stdlib.h>
    6 #include <iostream.h>
    7 #include <values.h>
    8 #include <complex>
    96#include "pexceptions.h"
    107#include "tmatrix.h"
    11 #include "objfio.h"
     8
     9
     10
    1211
    1312////////////////////////////////////////////////////////////////
    1413//**** Createur, Destructeur
    15 
    1614template <class T>
    1715TMatrix<T>::TMatrix()
    1816// Constructeur par defaut.
    19 : mNr(0), mNc(0), mNDBlock()
    20 {
    21 }
    22 
    23 template <class T>
    24 TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
     17  : TArray<T>()
     18{
     19}
     20
     21template <class T>
     22TMatrix<T>::TMatrix(uint_4 r,uint_4 c, short mm)
    2523// Construit une matrice de r lignes et c colonnes.
    26 : mNr(r), mNc(c), mNDBlock(r*c)
    27 {
    28 }
    29 
    30 template <class T>
    31 TMatrix<T>::TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br)
    32 // Construit une matrice de r lignes et c colonnes. On fournit
    33 // le tableau des valeurs et eventuellement un Bridge.
    34 : mNr(r), mNc(c), mNDBlock(r*c,values,br)
    35 {
     24  :  TArray<T>()
     25{
     26  if ( (r == 0) || (c == 0) )
     27    throw ParmError("TMatrix<T>::TMatrix(uint_4 r,uint_4 c) NRows or NCols = 0");
     28  ReSize(r, c, mm);
    3629}
    3730
    … …  
    3932TMatrix<T>::TMatrix(const TMatrix<T>& a)
    4033// Constructeur par copie (partage si "a" temporaire).
    41 : mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
    42 {
    43 }
    44 
    45 template <class T>
    46 TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
     34  : TArray<T>(a)
     35{
     36}
     37
     38template <class T>
     39TMatrix<T>::TMatrix(const TMatrix<T>& a, bool share)
    4740// Constructeur par copie avec possibilite de forcer le partage ou non.
    48 : mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
    49 {
     41: TArray<T>(a, share)
     42{
     43}
     44
     45
     46template <class T>
     47TMatrix<T>::TMatrix(const TArray<T>& a)
     48: TArray<T>(a)
     49{
     50  if (a.NbDimensions() != 2)
     51    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
     52}
     53
     54template <class T>
     55TMatrix<T>::TMatrix(const TArray<T>& a, bool share, short mm )
     56: TArray<T>(a, share)
     57{
     58  if (a.NbDimensions() != 2)
     59    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
     60  UpdateMemoryMapping(a, mm);
    5061}
    5162
    … …  
    5465// Destructeur
    5566{
    56 }
    57 
    58 ////////////////////////////////////////////////////////////////
    59 // Operations matricielles
    60 template <class T>
    61 TMatrix<T> TMatrix<T>::Transpose(void) const
     67
     68}
     69
     70template <class T>
     71TArray<T>& TMatrix<T>::Set(const TArray<T>& a)
     72{
     73  if (a.NbDimensions() != 2)
     74    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() != 2 ");
     75  return(TArray<T>::Set(a));
     76}
     77
     78template <class T>
     79void TMatrix<T>::ReSize(uint_4 r, uint_4 c, short mm)
     80{
     81  if(r==0||c==0)
     82    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
     83  uint_4 size[BASEARRAY_MAXNDIMS];
     84  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
     85  size[0] = r;  size[1] = c;
     86  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     87  UpdateMemoryMapping(mm, size[0], size[1]);
     88  TArray<T>::ReSize(2, size, 1);
     89  UpdateMemoryMapping(mm, size[0], size[1]);
     90}
     91
     92template <class T>
     93void TMatrix<T>::Realloc(uint_4 r,uint_4 c, short mm, bool force)
     94{
     95  if(r==0||c==0)
     96    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
     97  uint_4 size[BASEARRAY_MAXNDIMS];
     98  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
     99  UpdateMemoryMapping(mm, size[0], size[1]);
     100  TArray<T>::Realloc(2, size, 1, force);
     101  UpdateMemoryMapping(mm, size[0], size[1]);
     102}
     103
     104// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     105template <class T>
     106TMatrix<T> TMatrix<T>::operator () (Range rline, Range rcol) const
     107{
     108  uint_4 nx, ny;
     109  nx = ny = 0;
     110  if (GetMemoryMapping() == CMemoryMapping ) {
     111    TMatrix sm(SubArray(rcol, rline, 0, 0, 0),true,CMemoryMapping);
     112    sm.UpdateMemoryMapping(CMemoryMapping, nx, ny);
     113    sm.SetTemp(true);
     114    return(sm);
     115  }
     116  else {
     117    TMatrix sm(SubArray(rline, rcol, 0, 0, 0),true,CMemoryMapping);
     118    sm.UpdateMemoryMapping(FortranMemoryMapping, nx, ny);
     119    sm.SetTemp(true);
     120    return(sm);
     121  }
     122}
     123
     124////////////////////////////////////////////////////////////////
    62125// Transposition
    63 {
    64 TMatrix<T> a; a.Clone(*this); a.SetTemp(true);
    65 a.mNr = mNc; a.mNc = mNr;
    66 {for(uint_4 i=0; i<mNr; i++)
    67   for(uint_4 j=0; j<mNc; j++) {
    68     a(j,i) = (*this)(i,j);
    69 }}
    70 return a;
    71 }
     126template <class T>
     127TMatrix<T>& TMatrix<T>::Transpose()
     128{
     129  uint_4 rci = macoli_;
     130  macoli_ = marowi_;
     131  marowi_ = rci;
     132  return(*this);
     133}
     134
     135
     136template <class T>
     137TMatrix<T> TMatrix<T>::Transpose(short mm)
     138{
     139  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     140  TMatrix<T> tm(NCols(), NRows(), mm);
     141  for(uint_4 i=0; i<NRows(); i++)
     142    for(uint_4 j=0; j<NCols(); j++)
     143      tm(j,i) = (*this)(i,j);
     144  tm.SetTemp(true);
     145  return tm;
     146}
     147
     148// Rearrangement memoire
     149template <class T>
     150TMatrix<T> TMatrix<T>::Rearrange(short mm)
     151{
     152  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     153  TMatrix<T> tm(NRows(), NCols(), mm);
     154  for(uint_4 i=0; i<NRows(); i++)
     155    for(uint_4 j=0; j<NCols(); j++)
     156      tm(i,j) = (*this)(i,j);
     157  tm.SetTemp(true);
     158  return tm;
     159}
     160
     161template <class T>
     162TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
     163{
     164  if (ndim_ == 0) {
     165    uint_4 sz = imx.Size();
     166    if (sz < 1) sz = 1;
     167    ReSize(sz, sz);
     168  }
     169  T diag = (T)imx.Diag();
     170  if (NRows() != NCols())
     171    throw SzMismatchError("TMatrix::operator= (IdentityMatrix) NRows() != NCols()") ;
     172  for(uint_4 i=0; i<NRows(); i++) (*this)(i,i) = diag;
     173
     174  return (*this);
     175}
     176
     177
    72178
    73179////////////////////////////////////////////////////////////////
    74180//**** Impression
    75 
    76 template <class T>
    77 void TMatrix<T>::Print(ostream& os,int lp
    78                  ,uint_4 i0,uint_4 ni,uint_4 j0,uint_4 nj) const
    79 // Impression de la sous-matrice (i:i+ni-1,i:j+nj-1)
    80 {
    81 os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
    82 if(lp>0)
    83   {os<<"       this="<<this<<endl; mNDBlock.Print(0,0);}
    84 if(mNr==0 || mNc==0) return;
    85 if(i0>=mNr || j0>=mNc || ni==0 || nj==0) return;
    86 uint_4 i1 = i0+ni; if(i1>mNr) i1 = mNr;
    87 uint_4 j1 = j0+nj; if(j1>mNc) j1 = mNc;
    88 for(uint_4 i=i0;i<i1;i++) {
    89   for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
    90   cout<<endl;
    91 }
    92 }
    93 
    94 ////////////////////////////////////////////////////////////////
    95 //**** Surcharge de *= (INPLACE): TMatrix *= TMatrix;
    96 
    97 template <class T>
    98 TMatrix<T>& TMatrix<T>::operator *= (const TMatrix<T>& a)
    99 // A = A*B  -> A(n,m) = A(n,m)*B(m,m)
    100 {
    101 uint_4 ndata = mNr*mNc;
    102 if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc)
    103   throw(SzMismatchError("TMatrix::operator*=A size mismatch"));
    104 // A(i,j) = Sum(k) A(i,k)*B(k,j) ... il faut sauver la ligne "i" de A
    105 // Vecteur oi   : vecteur ou est sauve la ligne i de la matrice *this
    106 //                oi,oe = pointeur de debut et de fin du vecteur temporaire
    107 //                oij   = pointeur parcourant le vecteur oi
    108 // Matrice *this: i     = pointeur du debut de la ligne i
    109 //                ij    = pointeur parcourant la ligne i
    110 // Matrice a    : aj    = pointeur de debut de la colonne j
    111 //                aji   = pointeur parcourant la colonne j
    112 T* oi = new T[mNc]; T* oe = oi+mNc;
    113 for(T *i=Data(); i<Data()+ndata; i+=mNc) {
    114   {for(T *oij=oi,*ij=i; oij<oe;) *oij++ = *ij++;}
    115   {for(T *ij=i,*aj=const_cast<T *>(a.Data()); aj<a.Data()+a.mNc; ij++,aj++) {
    116     T sum = 0;
    117     for(T *oij=oi,*aji=aj; oij<oe; oij++,aji+=a.mNc) sum += *oij * *aji;
    118     *ij = sum;
    119   }}
    120 }
    121 delete [] oi;
    122 return *this;
    123 }
    124 
    125 ////////////////////////////////////////////////////////////////
    126 //**** Pour surcharge d'operateurs C = A (+,-,*,/) B
    127 
    128 template <class T> TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
    129 {
    130 if(mNr!=b.mNr || mNc!=b.mNc)
    131     throw(SzMismatchError("TMatrix operator C=A+B size mismatch\n"));
    132 TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
    133 result.mNDBlock = mNDBlock+b.mNDBlock;
    134 return result;
    135 }
    136 
    137 template <class T> TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
    138 {
    139 if(mNr!=b.mNr || mNc!=b.mNc)
    140     throw(SzMismatchError("TMatrix operator C=A-B size mismatch\n"));
    141 TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
    142 result.mNDBlock = mNDBlock-b.mNDBlock;
    143 return result;
    144 }
    145 
    146 template <class T> TMatrix<T> TMatrix<T>::Mul(const TMatrix<T>& b) const
    147 // C = A(this)*B : Cij = Aik Bkj (allocation forcee dans tous les cas)
    148 {
    149 if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
    150   throw(SzMismatchError("TMatrix operator C=A*B size mismatch\n"));
    151 TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
    152 T *ai,*aik,*bj,*bkj,*ri,*rij;
    153 for(ri=const_cast<T *>(r.Data()),ai=const_cast<T *>(Data());
    154     ri<r.Data()+r.mNr*r.mNc;ri+=r.mNc,ai+=mNc) {
    155   for(rij=ri,bj=const_cast<T *>(b.Data());rij<ri+r.mNc;rij++,bj++) {
    156     *rij = 0;
    157     for(aik=ai,bkj=bj;aik<ai+mNc;aik++,bkj+=b.mNc) *rij += *aik * *bkj;
    158   }
    159 }
    160 return r;
    161 }
    162 
    163 
    164 ///////////////////////////////////////////////////////////
    165 // --------------------------------------------------------
    166 //   Les objets delegues pour la gestion de persistance
    167 // --------------------------------------------------------
    168 ///////////////////////////////////////////////////////////
    169 
    170 template <class T>
    171 FIO_TMatrix<T>::FIO_TMatrix()
    172 {
    173 dobj=new TMatrix<T>;
    174 ownobj=true;
    175 }
    176 
    177 template <class T>
    178 FIO_TMatrix<T>::FIO_TMatrix(string const & filename)
    179 {
    180 dobj=new TMatrix<T>;
    181 ownobj=true;
    182 Read(filename);
    183 }
    184 
    185 template <class T>
    186 FIO_TMatrix<T>::FIO_TMatrix(const TMatrix<T> & obj)
    187 {
    188 dobj = new TMatrix<T>(obj);
    189 ownobj=true;
    190 }
    191 
    192 template <class T>
    193 FIO_TMatrix<T>::FIO_TMatrix(TMatrix<T> * obj)
    194 {
    195 dobj = obj;
    196 ownobj=false;
    197 }
    198 
    199 template <class T>
    200 FIO_TMatrix<T>::~FIO_TMatrix()
    201 {
    202 if (ownobj && dobj) delete dobj;
    203 }
    204 
    205 template <class T>
    206 AnyDataObj* FIO_TMatrix<T>::DataObj()
    207 {
    208 return(dobj);
    209 }
    210 
    211 template <class T>
    212 void FIO_TMatrix<T>::SetDataObj(AnyDataObj & o)
    213 {
    214 TMatrix<T> * po = dynamic_cast< TMatrix<T> * >(&o);
    215 if (po == NULL) return;
    216 if (ownobj && dobj) delete dobj;
    217 dobj = po; ownobj = false;
    218 }
    219 
    220 template <class T>
    221 void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
    222 {
    223 // On lit les 3 premiers uint_4
    224 //  0: Numero de version,  1 : NRows,  2 : NCol
    225 uint_4 itab[3];
    226 is.Get(itab,3);
    227 if (dobj == NULL) dobj = new TMatrix<T>(itab[1],itab[2]);
    228 else dobj->ReSize(itab[1],itab[2]);
    229 // On lit le NDataBlock
    230 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
    231 fio_nd.Read(is);
    232 }
    233 
    234 template <class T>
    235 void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
    236 {
    237 if (dobj == NULL)   return;
    238 //  On ecrit 3 uint_4 ....
    239 //  0: Numero de version,  1 : NRows,  2 : NCol
    240 uint_4 itab[3];
    241  itab[0] = 1;  // Numero de version a 1
    242 itab[1] = dobj->NRows();
    243 itab[2] = dobj->NCols();
    244 os.Put(itab,3);
    245 // On ecrit le NDataBlock
    246 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
    247 fio_nd.Write(os);
    248 }
    249 
    250 ////////////////////////////////////////////////////////////////
    251 // -------------------------------------------------------------
    252 //   La classe de gestion des lignes et colonnes d'une matrice
    253 // -------------------------------------------------------------
    254 ////////////////////////////////////////////////////////////////
    255 #include "tvector.h"
     181template <class T>
     182void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
     183{
     184  os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols()
     185     << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
     186  os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
     187  if (maxprt < 0)  maxprt = max_nprt_;
     188  int_4 npr = 0;
     189  Show(os, si);
     190  uint_4 kc,kr; 
     191  for(kr=0; kr<size_[marowi_]; kr++) {
     192    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) ) cout << "----- Ligne Line= " << kr << endl;
     193    for(kc=0; kc<size_[macoli_]; kc++) {
     194      if(kc > 0) os << ", "; 
     195      os << (*this)(kr, kc);   npr++;
     196      if (npr >= maxprt) {
     197        if (npr < totsize_)  os << "\n     .... " << endl; return;
     198      }
     199    }
     200    os << endl;
     201  }
     202  os << "\n" << endl;
     203}
     204
     205////////////////////////////////////////////////////////////////
     206//****  Multiplication matricielle  *****
     207////////////////////////////////////////////////////////////////
     208
     209template <class T>
     210TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
     211{
     212  if (NCols() != b.NRows())
     213    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
     214  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     215  TMatrix<T> rm(NRows(), b.NCols(), mm);
     216
     217  const T * pea;
     218  const T * peb;
     219  T sum;
     220  uint_4 r,c,k;
     221  uint_4 stepa = Step(ColsKA());
     222  uint_4 stepb = b.Step(RowsKA());
     223  // Calcul de C=rm = A*B   (A=*this)
     224  for(r=0; r<rm.NRows(); r++)      // Boucle sur les lignes de A
     225    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
     226      sum = 0;
     227      pea = &((*this)(r,0));       // 1er element de la ligne r de A
     228      peb = &(b(0,c));                // 1er element de la colonne c de B
     229      for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
     230      rm(r,c) = sum;
     231    }
     232
     233  rm.SetTemp(true);
     234  return rm;
     235}
    256236
    257237
    258238///////////////////////////////////////////////////////////////
    259239#ifdef __CXX_PRAGMA_TEMPLATES__
    260 #pragma define_template TMatrix<uint_1>
    261240#pragma define_template TMatrix<uint_2>
    262 #pragma define_template TMatrix<int_2>
    263241#pragma define_template TMatrix<int_4>
    264242#pragma define_template TMatrix<int_8>
    265 #pragma define_template TMatrix<uint_4>
    266 #pragma define_template TMatrix<uint_8>
    267243#pragma define_template TMatrix<r_4>
    268 #pragma define_template TMatrix<r_8>
     244#pragma define_template TMatrix<r_8> 
    269245#pragma define_template TMatrix< complex<r_4> >
    270246#pragma define_template TMatrix< complex<r_8> >
    271 // Instances des delegues FileIO (PPersist)
    272 #pragma define_template FIO_TMatrix<uint_1>
    273 #pragma define_template FIO_TMatrix<uint_2>
    274 #pragma define_template FIO_TMatrix<int_2>
    275 #pragma define_template FIO_TMatrix<int_4>
    276 #pragma define_template FIO_TMatrix<int_8>
    277 #pragma define_template FIO_TMatrix<uint_4>
    278 #pragma define_template FIO_TMatrix<uint_8>
    279 #pragma define_template FIO_TMatrix<r_8>
    280 #pragma define_template FIO_TMatrix<r_4>
    281 #pragma define_template FIO_TMatrix< complex<r_4> >
    282 #pragma define_template FIO_TMatrix< complex<r_8> >
    283247#endif
    284248
    285249#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    286 template class TMatrix<uint_1>;
    287250template class TMatrix<uint_2>;
    288 template class TMatrix<int_2>;
    289251template class TMatrix<int_4>;
    290252template class TMatrix<int_8>;
    291 template class TMatrix<uint_4>;
    292 template class TMatrix<uint_8>;
    293253template class TMatrix<r_4>;
    294254template class TMatrix<r_8>;
    295255template class TMatrix< complex<r_4> >;
    296256template class TMatrix< complex<r_8> >;
    297 // Instances des delegues FileIO (PPersist)
    298 template class FIO_TMatrix<uint_1>;
    299 template class FIO_TMatrix<uint_2>;
    300 template class FIO_TMatrix<int_2>;
    301 template class FIO_TMatrix<int_4>;
    302 template class FIO_TMatrix<int_8>;
    303 template class FIO_TMatrix<uint_4>;
    304 template class FIO_TMatrix<uint_8>;
    305 template class FIO_TMatrix<r_8>;
    306 template class FIO_TMatrix<r_4>;
    307 template class FIO_TMatrix< complex<r_4> >;
    308 template class FIO_TMatrix< complex<r_8> >;
    309257#endif

  • trunk/SophyaLib/TArray/tmatrix.h

    r772 r804  
    55
    66#include "machdefs.h"
    7 #include <stdio.h>
    8 #include <iostream.h>
    9 #include <complex>
    10 #include "ppersist.h"
    11 #include "anydataobj.h"
    12 #include "ndatablock.h"
     7#include "tarray.h"
    138
    149namespace SOPHYA {
    1510
    1611template <class T>
    17 class TMatrix : public AnyDataObj {
     12class TMatrix : public TArray<T> {
    1813public:
    19 
    2014  // Creation / destruction
    2115  TMatrix();
    22   TMatrix(uint_4 r,uint_4 c);
    23   TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
     16  TMatrix(uint_4 r,uint_4 c, short mm=AutoMemoryMapping);
    2417  TMatrix(const TMatrix<T>& a);
    25   TMatrix(const TMatrix<T>& a,bool share);
     18  TMatrix(const TMatrix<T>& a, bool share);
     19  TMatrix(const TArray<T>& a);
     20  TMatrix(const TArray<T>& a,  bool share, short mm=CMemoryMapping);
    2621  virtual ~TMatrix();
    2722
    28   // Temporaire?
    29   inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
    30   inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
     23  //  Pour verifiez la compatibilite de dimensions lors de l'affectation
     24  virtual TArray<T>& Set(const TArray<T>& a);
     25  inline  TMatrix<T>& operator = (const TMatrix<T>& a)
     26                       { Set(a);  return(*this); }
    3127
    32   // Gestion taille/Remplissage
    33   inline void Clone(const TMatrix<T>& a)  // Clone: copie des donnees de "a"
    34               {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
    35   inline void Reset(T v=0) {mNDBlock.Reset(v);}
    36   inline void ReSize(uint_4 r,uint_4 c)  // Reallocation de place
    37        {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
    38         mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
    39   inline void Realloc(uint_4 r,uint_4 c,bool force=false)
    40        {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
    41         mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}
     28  // Size - Changing the Size
     29  inline uint_4 NRows() const {return Size(marowi_); }
     30  inline uint_4 NCols() const {return Size(macoli_); }
     31  inline uint_4 NCol()  const {return Size(macoli_); } // back-compat Peida
    4232
    43   // Informations pointeur/data
    44   inline uint_4 NRows() const {return mNr;}
    45   inline uint_4 NCols() const {return mNc;}
    46   inline uint_4 NCol() const {return mNc;} // back-compat Peida
    47   inline T const& operator()(uint_4 r,uint_4 c) const
    48                             {return *(mNDBlock.Begin()+r*mNc+c);}
    49   inline T&       operator()(uint_4 r,uint_4 c)
    50                             {return *(mNDBlock.Begin()+r*mNc+c);}
    51   inline T const& operator[](uint_4 ip) const
    52                             {return *(mNDBlock.Begin()+ip);}
    53   inline T&       operator[](uint_4 ip)
    54                             {return *(mNDBlock.Begin()+ip);}
    55   inline       T* Data()       {return mNDBlock.Begin();}
    56   inline const T* Data() const {return mNDBlock.Begin();}
    57   inline       NDataBlock<T>& DataBlock()       {return mNDBlock;}
    58   inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}
     33  void ReSize(uint_4 r,uint_4 c, short mm=SameMemoryMapping);  // Reallocation de place
     34  void Realloc(uint_4 r,uint_4 c, short mm=SameMemoryMapping, bool force=false);
     35
     36  // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     37  TMatrix<T> operator () (Range rline, Range rcol) const ;
     38
     39  // Inline element acces methods
     40  inline T const& operator()(uint_4 r,uint_4 c) const;
     41  inline T&       operator()(uint_4 r,uint_4 c);
    5942
    6043  // Operations matricielles
    61   TMatrix<T> Transpose(void) const;
     44  TMatrix<T>& Transpose();
     45  //mm = SameMemoryMapping or CMemoryMapping or FortranMemoryMapping
     46  TMatrix<T>  Transpose(short mm);
     47  // Rearranging Matrix Elements
     48  TMatrix<T>  Rearrange(short mm);
    6249
    6350  // Operateur d'affectation
    64   // A = x (matrice diagonale x*Identite)
    65   inline TMatrix<T>& operator = (T x)
    66     {if(mNr!=mNc || mNr==0) throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
    67      for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++) (*this)(r,c)=(r==c)?x:(T)0;
    68      return *this;}
    69   // A = B : partage les donnees si "a" est temporaire, clone sinon.
    70   inline TMatrix<T>& operator = (const TMatrix<T>& a)
    71                 {if(this == &a) return *this; CloneOrShare(a); return *this;}
     51  // A = x (matrice diagonale Identite)
     52  virtual TMatrix<T>& SetIdentity(IdentityMatrix imx);
     53  inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }
    7254
    73   // Impression
    74   void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
    75   inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
    76               {Print(cout,lp,i0,ni,j0,nj);}
     55  // Operations diverses  avec une constante
     56  inline  TMatrix<T>&  operator = (T x)             { Set(x); return(*this); }
     57  inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
     58  inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
     59  inline  TMatrix<T>&  operator *= (T x)            { Mul(x); return(*this); }
     60  inline  TMatrix<T>&  operator /= (T x)            { Div(x); return(*this); }
    7761
    78   // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
    79   inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
    80   inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
    81   inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
    82   inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}
     62  //  operations avec matrices
     63  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { return AddElt(a); }
     64  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { return SubElt(a); }
    8365
    84   // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
    85   inline TMatrix<T>& operator += (const TMatrix<T>& a)
    86               {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
    87                   throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
    88                mNDBlock += a.mNDBlock; return *this;}
    89   inline TMatrix<T>& operator -= (const TMatrix<T>& a)
    90               {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
    91                   throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
    92                mNDBlock -= a.mNDBlock; return *this;}
    93   TMatrix<T>& operator *= (const TMatrix<T>& a);
     66  // Input-Output
     67  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
    9468
    95   // Pour surcharge d'operateurs C = A (+,-,*) B
    96   TMatrix<T> Add(const TMatrix<T>& b) const;
    97   TMatrix<T> Sub(const TMatrix<T>& b) const;
    98   TMatrix<T> Mul(const TMatrix<T>& b) const;
     69  // Produit matriciel
     70  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
     71  //  inline TMatrix<T>& operator *= (const TMatrix<T>& b)
     72
     73protected:
     74};
     75
     76//  ---- inline acces methods ------
     77template <class T>
     78inline T const& TMatrix<T>::operator()(uint_4 r, uint_4 c) const
     79{
     80#ifdef SO_BOUNDCHECKING
     81  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
     82  else   CheckBound(c, r, 0, 0, 0, 4);
     83#endif
     84  return ( *( mNDBlock.Begin()+ offset_+
     85              r*step_[marowi_] + c*step_[macoli_] ) );
     86}
     87
     88template <class T>
     89inline T & TMatrix<T>::operator()(uint_4 r, uint_4 c)
     90{
     91#ifdef SO_BOUNDCHECKING
     92  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
     93  else   CheckBound(c, r, 0, 0, 0, 4);
     94#endif
     95  return ( *( mNDBlock.Begin()+ offset_+
     96              r*step_[marowi_] + c*step_[macoli_] ) );
     97}
    9998
    10099
    101 protected:
    102   // partage les donnees si "a" temporaire, clone sinon.
    103   inline void CloneOrShare(const TMatrix<T>& a)
    104               {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
    105   // Share: partage les donnees de "a"
    106   inline void Share(const TMatrix<T>& a)
    107               {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
     100// Surcharge d'operateurs C = A * B
    108101
    109   uint_4 mNr,mNc;
    110   NDataBlock<T> mNDBlock;
    111 };
     102template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, const TMatrix<T>& b)
     103{ TMatrix<T> result(a); result.SetTemp(true); result.Multiply(b);  return result;}
    112104
    113 ////////////////////////////////////////////////////////////////
    114 // Impression
    115 
    116 template <class T>
    117 inline ostream& operator << (ostream& os, const TMatrix<T>& a)
    118                             {a.Print(os); return(os);}
    119 
    120 ////////////////////////////////////////////////////////////////
    121 // Surcharge d'operateurs A (+,-,*,/) (T) x
    122 
    123 template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
    124     {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
    125 
    126 template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
    127     {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
    128 
    129 template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
    130     {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}
    131 
    132 template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
    133     {TMatrix<T> result(a); result.SetTemp(true);
    134      result.DataBlock() = b-result.DataBlock(); return result;}
    135 
    136 template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
    137     {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
    138 
    139 template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
    140     {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
    141 
    142 template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
    143     {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}
    144 
    145 ////////////////////////////////////////////////////////////////
    146 // Surcharge d'operateurs C = A (+,-,*,/) B
    147 
    148 template <class T>
    149 inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
    150                   {return a.Add(b);}
    151 
    152 template <class T>
    153 inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
    154                   {return a.Sub(b);}
    155 
    156 template <class T>
    157 inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
    158                   {return a.Mul(b);}
    159 
    160 ////////////////////////////////////////////////////////////////
    161 // Typedef pour simplifier et compatibilite Peida
    162105typedef TMatrix<r_8> Matrix;
    163 
    164 /////////////////////////////////////////////////////////////////////////
    165 // Classe pour la gestion de persistance
    166 template <class T>
    167 class FIO_TMatrix : public  PPersist  {
    168 public:
    169   FIO_TMatrix();
    170   FIO_TMatrix(string const & filename);
    171   FIO_TMatrix(const TMatrix<T> & obj);
    172   FIO_TMatrix(TMatrix<T> * obj);
    173   virtual ~FIO_TMatrix();
    174   virtual AnyDataObj* DataObj();
    175   virtual void        SetDataObj(AnyDataObj & o);
    176   inline operator TMatrix<T>() { return(*dobj); }
    177 protected :
    178   virtual void ReadSelf(PInPersist&);           
    179   virtual void WriteSelf(POutPersist&) const; 
    180   TMatrix<T> * dobj;
    181   bool ownobj;
    182 };
    183 
    184 template <class T>
    185 inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
    186 { FIO_TMatrix<T> fio(&obj);  fio.Write(os);  return(os); }
    187 template <class T>
    188 inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
    189 { FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }
    190 
    191106
    192107} // Fin du namespace

  • trunk/SophyaLib/TArray/tvector.cc

    r762 r804  
    1 // $Id: tvector.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
     1// $Id: tvector.cc,v 1.2 2000-04-03 17:35:59 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    4 #include <stdio.h>
    54#include <stdlib.h>
    6 #include <iostream.h>
    7 #include <complex>
    85#include "pexceptions.h"
    96#include "tvector.h"
    10 #include "objfio.h"
     7
     8////////////////////////////////////////////////////////////////
     9//**** Createur, Destructeur
    1110
    1211template <class T>
    13 TVector<T>::TVector(uint_4 n)
    14 // Constructeur d'un vecteur de "n" elements
    15 : TMatrix<T>(n,1)
     12TVector<T>::TVector()
     13  : TMatrix<T>()
    1614{
    1715}
    1816
    1917template <class T>
    20 TVector<T>::TVector(uint_4 n, T* values,Bridge* br)
    21 // Construit un vecteur de n elements. On fournit
    22 // le tableau des valeurs et eventuellement un Bridge.
    23 : TMatrix<T>(n,1,values,br)
     18TVector<T>::TVector(uint_4 n, short lcv, short mm)
     19// Constructeur
     20  : TMatrix<T>((lcv == ColumnVector) ? n : 1,(lcv == ColumnVector) ? 1 : 0, mm)
     21{
     22  if (lcv == ColumnVector) veceli_ = marowi_;
     23  else veceli_ = macoli_;
     24}
     25
     26
     27template <class T>
     28TVector<T>::TVector(const TVector<T>& a)
     29// Constructeur par copie (partage si "a" temporaire).
     30  : TMatrix<T>(a)
    2431{
    2532}
    2633
    2734template <class T>
    28 TVector<T>::TVector(const TVector<T>& v)
    29 // Constructeur par copie (partage si "v" temporaire).
    30 : TMatrix<T>(v)
     35TVector<T>::TVector(const TVector<T>& a, bool share)
     36// Constructeur par copie avec possibilite de forcer le partage ou non.
     37: TMatrix<T>(a, share)
    3138{
    3239}
    3340
    3441template <class T>
    35 TVector<T>::TVector(const TVector<T>& v,bool share)
    36 // Constructeur par copie avec possibilite de forcer le partage ou non.
    37 : TMatrix<T>(v,share)
     42TVector<T>::TVector(const TArray<T>& a)
     43: TMatrix<T>(a)
    3844{
    3945}
    4046
     47
    4148template <class T>
    42 TVector<T>::TVector(const TMatrix<T>& a)
    43 // Constructeur a partir d'une matrice "a" de dimension (n,1)
    44 : TMatrix<T>(a)
     49TVector<T>::TVector(const TArray<T>& a, bool share, short lcv, short mm )
     50: TMatrix<T>(a, share, mm)
    4551{
    46 if(a.NCols() != 1)
    47   throw(SzMismatchError("TVector(const TMatrix<T>& a) size mismatch, ncol!=1"));
     52  if (lcv == SameTypeVector) veceli_ = a.VectKA();
     53  else {
     54    if (lcv == ColumnVector) veceli_ = marowi_;
     55    else veceli_ = macoli_;
     56  }
     57}
     58
     59template <class T>
     60TVector<T>::~TVector()
     61// Destructeur
     62{
     63
     64}
     65
     66template <class T>
     67void TVector<T>::ReSize(uint_4 n, short lcv)
     68{
     69  if( n == 0 )
     70    throw(SzMismatchError("TVector::ReSize() n = 0 "));
     71  uint_4 r,c;
     72  if (lcv == SameTypeVector)  lcv = GetVectorType();
     73  if (lcv == ColumnVector) { r = n;  c = 1; }
     74  else { c = n; r = 1; }
     75  TMatrix<T>::ReSize(r,c);
     76}
     77
     78template <class T>
     79void TVector<T>::Realloc(uint_4 n, short lcv, bool force)
     80{
     81  if( n == 0 )
     82    throw(SzMismatchError("TVector::Realloc() n = 0 "));
     83  uint_4 r,c;
     84  if (lcv == SameTypeVector)  lcv = GetVectorType();
     85  if (lcv == ColumnVector) { r = n;  c = 1; }
     86  else { c = n; r = 1; }
     87  TMatrix<T>::Realloc(r,c,SameMemoryMapping,force);
     88}
     89
     90// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     91template <class T>
     92TVector<T> TVector<T>::operator () (Range relt) const
     93{
     94  Range rr, cr;
     95  if (GetVectorType() == ColumnVector )  rr = relt;
     96  else cr = relt;
     97  TMatrix<T> const & mtx = (*this);
     98  TVector sv( mtx(rr, cr) , true, GetVectorType(), GetMemoryMapping() );
     99  sv.SetTemp(true);
     100  return(sv); 
     101}
     102
     103template <class T>
     104T TVector<T>::Norm2() const
     105{
     106  T ret = 0;
     107  for(uint_8 k=0; k<Size(); k++)   ret += (*this)(k)*(*this)(k);
     108  return ret;
    48109}
    49110
    50111
    51 ///////////////////////////////////////////////////////////
    52 // --------------------------------------------------------
    53 //   Les objets delegues pour la gestion de persistance
    54 // --------------------------------------------------------
    55 ///////////////////////////////////////////////////////////
    56 
    57 template <class T>
    58 FIO_TVector<T>::FIO_TVector()
    59 {
    60 dobj=new TVector<T>;
    61 ownobj=true;
    62 }
    63 
    64 template <class T>
    65 FIO_TVector<T>::FIO_TVector(string const & filename)
    66 {
    67 dobj=new TVector<T>;
    68 ownobj=true;
    69 Read(filename);
    70 }
    71 
    72 template <class T>
    73 FIO_TVector<T>::FIO_TVector(const TVector<T> & obj)
    74 {
    75 dobj = new TVector<T>(obj);
    76 ownobj=true;
    77 }
    78 
    79 template <class T>
    80 FIO_TVector<T>::FIO_TVector(TVector<T> * obj)
    81 {
    82 dobj = obj;
    83 ownobj=false;
    84 }
    85 
    86 template <class T>
    87 FIO_TVector<T>::~FIO_TVector()
    88 {
    89 if (ownobj && dobj) delete dobj;
    90 }
    91 
    92 template <class T>
    93 AnyDataObj* FIO_TVector<T>::DataObj()
    94 {
    95 return(dobj);
    96 }
    97 
    98 template <class T>
    99 void FIO_TVector<T>::SetDataObj(AnyDataObj & o)
    100 {
    101 TVector<T> * po = dynamic_cast< TVector<T> * >(&o);
    102 if (po == NULL) return;
    103 if (ownobj && dobj) delete dobj;
    104 dobj = po; ownobj = false;
    105 }
    106 
    107 template <class T>
    108 void FIO_TVector<T>::ReadSelf(PInPersist& is)
    109 {
    110 // On lit les 3 premiers uint_4
    111 //  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
    112 uint_4 itab[3];
    113 is.Get(itab,3);
    114 if (dobj == NULL) dobj = new TVector<T>(itab[1]);
    115 else dobj->ReSize(itab[1]);
    116 // On lit le NDataBlock
    117 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
    118 fio_nd.Read(is);
    119 }
    120 
    121 template <class T>
    122 void FIO_TVector<T>::WriteSelf(POutPersist& os) const
    123 {
    124 if (dobj == NULL)   return;
    125 //  On ecrit 3 uint_4 ....
    126 //  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
    127 uint_4 itab[3];
    128  itab[0] = 1;  // Numero de version a 1
    129 itab[1] = dobj->NElts();
    130 itab[2] = 1;
    131 os.Put(itab,3);
    132 // On ecrit le NDataBlock
    133 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
    134 fio_nd.Write(os);
    135 }
    136 
    137112///////////////////////////////////////////////////////////////
    138113#ifdef __CXX_PRAGMA_TEMPLATES__
    139 #pragma define_template TVector<uint_1>
    140114#pragma define_template TVector<uint_2>
    141 #pragma define_template TVector<int_2>
    142115#pragma define_template TVector<int_4>
    143116#pragma define_template TVector<int_8>
    144 #pragma define_template TVector<uint_4>
    145 #pragma define_template TVector<uint_8>
    146117#pragma define_template TVector<r_4>
    147118#pragma define_template TVector<r_8>
    148119#pragma define_template TVector< complex<r_4> >
    149120#pragma define_template TVector< complex<r_8> >
    150 // Instances des delegues FileIO (PPersist)
    151 #pragma define_template FIO_TVector<uint_1>
    152 #pragma define_template FIO_TVector<uint_2>
    153 #pragma define_template FIO_TVector<int_2>
    154 #pragma define_template FIO_TVector<int_4>
    155 #pragma define_template FIO_TVector<int_8>
    156 #pragma define_template FIO_TVector<uint_4>
    157 #pragma define_template FIO_TVector<uint_8>
    158 #pragma define_template FIO_TVector<r_8>
    159 #pragma define_template FIO_TVector<r_4>
    160 #pragma define_template FIO_TVector< complex<r_4> >
    161 #pragma define_template FIO_TVector< complex<r_8> >
    162121#endif
    163122
    164123#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    165 template class TVector<uint_1>;
    166124template class TVector<uint_2>;
    167 template class TVector<int_2>;
    168125template class TVector<int_4>;
    169126template class TVector<int_8>;
    170 template class TVector<uint_4>;
    171 template class TVector<uint_8>;
    172127template class TVector<r_4>;
    173128template class TVector<r_8>;
    174129template class TVector< complex<r_4> >;
    175130template class TVector< complex<r_8> >;
    176 // Instances des delegues FileIO (PPersist)
    177 template class FIO_TVector<uint_1>;
    178 template class FIO_TVector<uint_2>;
    179 template class FIO_TVector<int_2>;
    180 template class FIO_TVector<int_4>;
    181 template class FIO_TVector<int_8>;
    182 template class FIO_TVector<uint_4>;
    183 template class FIO_TVector<uint_8>;
    184 template class FIO_TVector<r_8>;
    185 template class FIO_TVector<r_4>;
    186 template class FIO_TVector< complex<r_4> >;
    187 template class FIO_TVector< complex<r_8> >;
    188131#endif
     132

  • trunk/SophyaLib/TArray/tvector.h

    r772 r804  
    1111class TVector : public TMatrix<T> {
    1212public:
     13  enum {SameTypeVector=0, ColumnVector=1, LineVector=2};
     14  // Creation / destruction
     15  TVector();
     16  TVector(uint_4 n, short lcv=ColumnVector, short mm=AutoMemoryMapping);
     17  TVector(const TVector<T>& v);
     18  TVector(const TVector<T>& v, bool share);
     19  TVector(const TArray<T>& a);
     20  TVector(const TArray<T>& a,  bool share, short lcv=ColumnVector, short mm=CMemoryMapping);
    1321
    14   // Creation / destruction
    15   TVector(uint_4 n=1);
    16   TVector(uint_4 n, T* values,Bridge* br=NULL);
    17   TVector(const TVector<T>& v);
    18   TVector(const TVector<T>& v,bool share);
    19   TVector(const TMatrix<T>& a);
     22  virtual ~TVector();
     23
     24  inline  TVector<T>& operator = (const TVector<T>& a)
     25                       { Set(a);  return(*this); }
     26
     27  // Vector type   Line or Column vector
     28  inline short GetVectorType() const
     29               { return((marowi_ == veceli_) ? ColumnVector : LineVector); }
    2030
    2131  // Gestion taille/Remplissage
    22   inline void ReSize(uint_4 n) {TMatrix<T>::ReSize(n,1);} // Reallocation de place
    23   inline void Realloc(uint_4 n,bool force=false) {TMatrix<T>::Realloc(n,1,force);}
     32  void ReSize(uint_4 n, short lcv=SameTypeVector );
     33  void Realloc(uint_4 n, short lcv=SameTypeVector, bool force=false);
     34
     35  // Sub-Vector extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
     36  TVector<T> operator () (Range relt) const ;
    2437
    2538  // Informations pointeur/data
    26   inline uint_4 NElts() const {return NRows();}
     39  inline uint_4 NElts() const {return Size(); }
    2740 
    28   // Acces aux elements
    29   inline T& operator()(uint_4 n) {return (*this)[n];}
    30   inline T const& operator()(uint_4 n) const {return (*this)[n];}
    31   inline T& Element(uint_4 n) {return (*this)[n];}
    32   inline T const& Element(uint_4 n) const {return (*this)[n];}
     41  // Inline element acces methods
     42  inline T const& operator()(uint_4 n) const;
     43  inline T&       operator()(uint_4 n);
    3344
    34   // Operateur d'affectation
    35   inline TVector& operator = (const TVector& v)
    36                   {TMatrix<T>::operator =(v); return *this;}
    37   inline TVector& operator = (T x)
    38                   {for(uint_4 i=0;i<NRows();i++) (*this)(i)=x; return *this;}
     45  // Operations diverses  avec une constante
     46  inline  TVector<T>&  operator = (T x)             { Set(x); return(*this); }
     47  inline  TVector<T>&  operator += (T x)            { Add(x); return(*this); }
     48  inline  TVector<T>&  operator -= (T x)            { Sub(x); return(*this); }
     49  inline  TVector<T>&  operator *= (T x)            { Mul(x); return(*this); }
     50  inline  TVector<T>&  operator /= (T x)            { Div(x); return(*this); }
    3951
     52  //  operations avec matrices
     53  inline  TVector<T>&  operator += (const TVector<T>& a)  { AddElt(a); return(*this); }
     54  inline  TVector<T>&  operator -= (const TVector<T>& a)  { SubElt(a); return(*this); }
    4055
     56  // Norme(^2)
     57  T Norm2() const ;
    4158};
    4259
    43 // produit scalaire, matrice*vecteur
     60//  ---- inline acces methods ------
    4461template <class T>
    45 inline T operator* (const TVector<T>& v1, const TVector<T>& v2)
    46   {if(v1.NRows() != v2.NRows())
    47      throw(SzMismatchError("TVector::operator*(TVector& v1,TVector v2) size mismatch"));
    48    T *p = const_cast<T *>(v1.Data()), *pEnd = p+v1.NElts(),
    49      *q = const_cast<T *>(v2.Data()), r = 0;
    50    while (p<pEnd) r += *p++ * *q++;
    51    return r;}
     62inline T const& TVector<T>::operator()(uint_4 n) const
     63{
     64#ifdef SO_BOUNDCHECKING
     65  if (veceli__ == 0)   CheckBound(n, 0, 0, 0, 0, 4);
     66  else   CheckBound(0, n, 0, 0, 0, 4);
     67#endif
     68  return ( *( mNDBlock.Begin()+ offset_ + n*step_[veceli_] ) );
     69}
    5270
    5371template <class T>
    54 inline TVector<T> operator* (const TMatrix<T>& a, const TVector<T>& b)
    55    {return TVector<T>(a * ((TMatrix<T> const&)(b)));}
    56 
     72inline T & TVector<T>::operator()(uint_4 n)
     73{
     74#ifdef SO_BOUNDCHECKING
     75  if (veceli__ == 0)   CheckBound(n, 0, 0, 0, 0, 4);
     76  else   CheckBound(0, n, 0, 0, 0, 4);
     77#endif
     78  return ( *( mNDBlock.Begin()+ offset_ + n*step_[veceli_] ) );
     79}
    5780
    5881// Typedef pour simplifier et compatibilite Peida
    5982typedef TVector<r_8> Vector;
    6083
    61 /////////////////////////////////////////////////////////////////////////
    62 // Classe pour la gestion de persistance
    63 template <class T>
    64 class FIO_TVector : public  PPersist  {
    65 public:
    66   FIO_TVector();
    67   FIO_TVector(string const & filename);
    68   FIO_TVector(const TVector<T> & obj);
    69   FIO_TVector(TVector<T> * obj);
    70   virtual ~FIO_TVector();
    71   virtual AnyDataObj* DataObj();
    72   virtual void        SetDataObj(AnyDataObj & o);
    73   inline operator TVector<T>() { return(*dobj); }
    74 protected :
    75   virtual void ReadSelf(PInPersist&);           
    76   virtual void WriteSelf(POutPersist&) const; 
    77   TVector<T> * dobj;
    78   bool ownobj;
    79 };
    80 
    81 template <class T>
    82 inline POutPersist& operator << (POutPersist& os, TVector<T> & obj)
    83 { FIO_TVector<T> fio(&obj);  fio.Write(os);  return(os); }
    84 template <class T>
    85 inline PInPersist& operator >> (PInPersist& is, TVector<T> & obj)
    86 { FIO_TVector<T> fio(&obj);  fio.Read(is);  return(is); }
    87 
    8884} // Fin du namespace
    8985

  • trunk/SophyaLib/TArray/utilarr.cc

    r785 r804  
    2727}
    2828
     29/*
    2930Range & Range::operator = (uint_4 start)
    3031{
    … …  
    3435  return (*this);
    3536}
     37*/
     38
     39
     40IdentityMatrix::IdentityMatrix(double diag, uint_4 n)
     41{
     42  size_ = n;
     43  diag_ = diag;
     44}

  • trunk/SophyaLib/TArray/utilarr.h

    r785 r804  
    1818class Sequence {
    1919public:
    20   Sequence (double start=0., double step=1., Arr_DoubleFunctionOfX f=NULL);
     20  explicit Sequence (double start=0., double step=1., Arr_DoubleFunctionOfX f=NULL);
    2121  inline double & Start() { return start_; }
    2222  inline double & Step() { return step_; }
    … …  
    2929class Range {
    3030public:
    31   Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
     31  explicit Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
    3232  inline uint_4 & Start()  {  return start_; }
    3333  inline uint_4 & Size()  {  return size_; }
    3434  inline uint_4 & Step()  {  return step_; }
    35   Range & operator = (uint_4 start);
    3635protected:
    3736  uint_4 start_, size_, step_;
     37};
     38
     39class IdentityMatrix {
     40public:
     41  explicit IdentityMatrix(double diag=1., uint_4 n=0);
     42  inline uint_4 Size() { return size_; }
     43  inline double Diag() { return diag_; }
     44protected:
     45  uint_4 size_;
     46  double diag_;
    3847};
    3948
Note: See TracChangeset for help on using the changeset viewer.