Changeset 813 in Sophya


Ignore:
Timestamp:
Apr 5, 2000, 5:44:19 PM (25 years ago)
Author:
ansari
Message:

Correction bug/amelioarions TArray,TMatrix,TVector - Reza 5/4/2000

Location:
trunk/SophyaLib/TArray
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/basarr.cc

    r804 r813  
    1414                                 "Stride(int )", "ElemCheckBound()", "operator()" };
    1515uint_4 BaseArray::max_nprt_ = 50;
     16uint_4 BaseArray::prt_lev_ = 0;
    1617short  BaseArray::default_memory_mapping = CMemoryMapping;
    17 
    18 // Methodes statiques globales
    19 void BaseArray::SetMaxPrint(uint_4 nprt)
     18short  BaseArray::default_vector_type = ColumnVector;
     19uint_8 BaseArray::openmp_size_threshold = 200000;
     20
     21// ------ Methodes statiques globales --------
     22
     23//  1/ Gestion des impressions
     24void BaseArray::SetMaxPrint(uint_4 nprt, uint_4 lev)
    2025{
    2126  max_nprt_ = nprt;
     27  prt_lev_ = (lev < 3) ? lev : 3;
     28}
     29
     30void BaseArray::SetOpenMPSizeThreshold(uint_8 thr)
     31{
     32  openmp_size_threshold = thr;
    2233}
    2334
    … …  
    3243short BaseArray::SetDefaultMemoryMapping(short mm)
    3344{
    34   default_memory_mapping = ( (mm == CMemoryMapping) ? CMemoryMapping : FortranMemoryMapping) ;
     45  default_memory_mapping = (mm != CMemoryMapping) ? FortranMemoryMapping : CMemoryMapping;
    3546  return default_memory_mapping;
    3647}
    3748
    38 
     49short BaseArray::SetDefaultVectorType(short vt)
     50{
     51  default_vector_type = (vt != ColumnVector) ? RowVector : ColumnVector ;
     52  return default_vector_type;
     53}
     54
     55// Memory organisation
    3956short BaseArray::SelectMemoryMapping(short mm)
    4057// si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
    … …  
    4461  else return (default_memory_mapping);
    4562}
    46 
    47 void BaseArray::UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny)
    48 {
     63short BaseArray::SelectVectorType(short vt)
     64{
     65  if ((vt == ColumnVector) || (vt == RowVector))  return(vt);
     66  else return(default_vector_type);
     67}
     68
     69void BaseArray::UpdateMemoryMapping(short mm)
     70{
     71  short vt = default_vector_type;
    4972  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
    5073  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   } 
     74    marowi_  = 1;  macoli_ = 0;
     75  }
     76  else {
     77    marowi_ = 0;  macoli_ = 1;
     78  }
     79
     80  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
     81    // Choix automatique Vecteur ligne ou colonne
     82    if ( size_[macoli_] == 1)  veceli_ = marowi_;
     83    else veceli_ = macoli_;
     84  }
     85  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
     86  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize 
    5887}
    5988
    6089void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
    6190{
    62   if (mm == SameMemoryMapping)
     91  short vt = default_vector_type;
     92  if (mm == SameMemoryMapping) {
    6393    mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
     94    vt = (a.marowi_ == a.veceli_) ? ColumnVector : RowVector;
     95  }
     96  else if (mm == AutoMemoryMapping)   mm = default_memory_mapping;
     97
    6498  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
    6599  if (mm == CMemoryMapping) {
    66     marowi_ = veceli_ = 1;  macoli_ = 0;
    67   }
    68   else {
    69     marowi_ = veceli_ = 0;  macoli_ = 1;
    70   } 
    71 }
    72 
     100    marowi_  = 1;  macoli_ = 0;
     101  }
     102  else {
     103    marowi_ = 0;  macoli_ = 1;
     104  }
     105  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
     106    // Choix automatique Vecteur ligne ou colonne
     107    if ( size_[macoli_] == 1)  veceli_ = marowi_;
     108    else veceli_ = marowi_;
     109  }
     110  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
     111  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize 
     112}
     113
     114void BaseArray::SetMemoryMapping(short mm)
     115{
     116  if (mm == SameMemoryMapping) return;
     117  if (mm == AutoMemoryMapping)  mm = default_memory_mapping;
     118  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) return;
     119  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
     120  if (mm == CMemoryMapping) {
     121    marowi_ =  1;  macoli_ = 0;
     122  }
     123  else {
     124    marowi_ =  0;  macoli_ = 1;
     125  }
     126  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
     127    // Choix automatique Vecteur ligne ou colonne
     128    if ( size_[macoli_] == 1)  veceli_ = marowi_;
     129    else veceli_ = macoli_;
     130  }
     131  else   veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
     132  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize 
     133}
     134
     135void BaseArray::SetVectorType(short vt)
     136{
     137  if (vt == SameVectorType) return;
     138  if (vt == AutoVectorType) vt =  default_vector_type;
     139  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
     140    // Choix automatique Vecteur ligne ou colonne
     141    if ( size_[macoli_] == 1)  veceli_ = marowi_;
     142    else veceli_ = macoli_;
     143  }
     144  else  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
     145  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize 
     146}
    73147
    74148// -------------------------------------------------------
    … …  
    87161  moystep_ = 0;
    88162  offset_ = 0;
    89   // Default for matrices : Lines = Y , Columns = X, Default Vectore= ColumnVector
    90   macoli_ = 0;
    91   veceli_ = marowi_ = 1;
     163  // Default for matrices : Memory organisation and Vector type
     164  if (default_memory_mapping == CMemoryMapping) {
     165    marowi_ =  1;  macoli_ = 0;
     166  }
     167  else {
     168    marowi_ =  0;  macoli_ = 1;
     169  }
     170  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
     171  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize 
    92172}
    93173
    … …  
    103183  for(int k=0; k<ndim_; k++)
    104184    if (size_[k] != a.size_[k])  return(false);
    105   if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
    106        (veceli_ != a.veceli_) )  return(false);
     185  //  $CHECK$   Reza doit-on verifier ca
     186  if (ck_memo_vt_ && a.ck_memo_vt_)
     187    if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
     188         (veceli_ != a.veceli_) )  return(false);
    107189  return(true);
    108190}
    … …  
    166248
    167249//  Acces lineaire aux elements ....  Calcul d'offset
     250// --------------------------------------------------
     251// Position de l'element 0 du vecteur i selon l'axe ka
     252// --------------------------------------------------
     253uint_8 BaseArray::Offset(uint_4 ka, uint_8 i) const
     254{
     255
     256  if ( (ndim_ < 1) || (i == 0) )  return(offset_);
     257  //#ifdef SO_BOUNDCHECKING
     258  if (ka >= ndim_)
     259    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Axe KA Error");
     260  if ( i*size_[ka] >= totsize_ ) 
     261    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Index Error");
     262  //#endif
     263  uint_4 idx[BASEARRAY_MAXNDIMS];
     264  int k;
     265  uint_8 rest = i;
     266  idx[ka] = 0;
     267  for(k=0; k<ndim_; k++) {
     268    if (k == ka) continue;
     269    idx[k] = rest%size_[k];   rest /= size_[k];
     270  }
     271  uint_8 off = offset_;
     272  for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
     273  return (off);
     274}
    168275
    169276uint_8 BaseArray::Offset(uint_8 ip) const
    170277{
    171 if (ip == 0)  return(offset_);
    172 uint_4 idx[BASEARRAY_MAXNDIMS];
    173 int k;
    174 uint_8 rest = ip;
    175 for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
    176 uint_8 off = offset_;
    177 for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
    178 return (off);
     278  if ( (ndim_ < 1) || (ip == 0) )  return(offset_);
     279  //#ifdef SO_BOUNDCHECKING
     280  if (ip >= totsize_)
     281    throw RangeCheckError("BaseArray::Offset(uint_8 ip) Out of range index ip");
     282  //#endif
     283
     284  uint_4 idx[BASEARRAY_MAXNDIMS];
     285  int k;
     286  uint_8 rest = ip;
     287  for(k=0; k<ndim_; k++) {
     288    idx[k] = rest%size_[k];   rest /= size_[k];
     289  }
     290  //#ifdef SO_BOUNDCHECKING
     291  if (rest != 0)
     292    throw PError("BaseArray::Offset(uint_8 ip) GUG !!! rest != 0");
     293  //#endif
     294//   if (rest != 0) cerr << " BUG ---- BaseArray::Offset( " << ip << " )" << rest << endl;
     295//   cerr << " DBG-Offset( " << ip << ")" ;
     296//   for(k=0; k<ndim_; k++) cerr << idx[k] << "," ;
     297//   cerr << " ZZZZ " << endl;
     298  uint_8 off = offset_;
     299  for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
     300  return (off);
    179301}
    180302
    … …  
    186308void BaseArray::Show(ostream& os, bool si) const
    187309{
    188   os << " TArray< " << DataType() << " >  NDim= " << ndim_
    189      << "(" << typeid(*this).name() << " )" << endl;
    190   os << "TotSize=" << totsize_ << " Size(X*Y*...)= " ;
     310  os << "\n--- " << InfoString() ;
     311  os << " ND=" << ndim_ << " SizeX*Y*...= " ;
    191312  for(int k=0; k<ndim_; k++) {
    192313    os << size_[k];
    193     if (k<ndim_-1)  os << " x ";
    194   }
    195   os << endl;
    196   os <<  " Step(X Y ...)= " ;
    197   for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
    198   os << " Offset= " << offset_  << "\n " << endl;
    199   if (si && (mInfo != NULL)) os << (*mInfo) << endl;
     314    if (k<ndim_-1)  os << "x";
     315  }
     316  os << " ---" << endl;
     317  if (prt_lev_ > 0) {
     318    os <<  " TotSize= " << totsize_ << " Step(X Y ...)="  ;
     319    for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
     320    os << " Offset= " << offset_  << endl;
     321  }
     322  if (prt_lev_ > 1) {
     323    os << " MemoryMapping=" << GetMemoryMapping() << " VecType= " << GetVectorType()
     324       << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA()
     325       << " VectKA=" << VectKA() << endl;
     326  }
     327  if (!si && (prt_lev_ < 2)) return;
     328  if (mInfo != NULL) os << (*mInfo) << endl;
    200329 
    201330}
    202331
     332string BaseArray::InfoString() const
     333{
     334  string rs = "BaseArray Type= ";
     335  rs +=  typeid(*this).name() ;
     336  return rs;
     337}
    203338
    204339DVList& BaseArray::Info()
    … …  
    238373  }
    239374  offset_ = offset;
    240   // Default for matrices : Lines = Y , Columns = X
    241   macoli_ = 0;
    242   marowi_ = veceli_ = 1;
     375  // Default for matrices : Memory organisation and Vector type
     376  if (default_memory_mapping == CMemoryMapping) {
     377    marowi_ =  1;  macoli_ = 0;
     378  }
     379  else {
     380    marowi_ =  0;  macoli_ = 1;
     381  }
     382  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
     383  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize 
    243384  // Update OK
    244385  ndim_ = ndim;
    … …  
    280421  minstep_ = minstep;
    281422  offset_ = offset;
    282   // Default for matrices : Lines = Y , Columns = X
    283   macoli_ = 0;
    284   marowi_ = veceli_ = 1;
     423  // Default for matrices : Memory organisation and Vector type
     424  if (default_memory_mapping == CMemoryMapping) {
     425    marowi_ =  1;  macoli_ = 0;
     426  }
     427  else {
     428    marowi_ =  0;  macoli_ = 1;
     429  }
     430  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
     431  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize 
    285432  // Update OK
    286433  ndim_ = ndim;
    … …  
    323470  marowi_ = a.marowi_;
    324471  veceli_ = a.veceli_;
     472  ck_memo_vt_ = a.ck_memo_vt_;
    325473  // Update OK
    326474  ndim_ = a.ndim_;

  • trunk/SophyaLib/TArray/basarr.h

    r804 r813  
    2222public:
    2323  //  To define Array / Matrix memory mapping
    24   enum { AutoMemoryMapping = -1, SameMemoryMapping = 0,
     24  enum MemoryMapping { AutoMemoryMapping = -1, SameMemoryMapping = 0,
    2525         CMemoryMapping = 1, FortranMemoryMapping = 2 };
     26  //  Vector type
     27  enum VectorType {AutoVectorType = -1, SameVectorType = 0,
     28                   ColumnVector = 1, RowVector = 2};
    2629
    27   //    Max Nb of printed elements
    28   static void    SetMaxPrint(uint_4 nprt=50);   
    29   //   Memory organisation (for matrices)
     30  //    Size threshold for parallel routine call
     31  static void    SetOpenMPSizeThreshold(uint_8 thr=200000);
     32  static inline uint_8 GetOpenMPSizeThreshold() { return openmp_size_threshold; }
     33  //    Max Nb of printed elements and print level
     34  static void    SetMaxPrint(uint_4 nprt=50, uint_4 lev=0);   
     35  static inline uint_4 GetMaxPrint() { return max_nprt_; }
     36  static inline uint_4 GetPrintLevel() { return prt_lev_; }
     37
     38  //   Memory organisation (for matrices) and vector type
    3039  static short   SetDefaultMemoryMapping(short mm=CMemoryMapping);
    3140  static inline short GetDefaultMemoryMapping() { return default_memory_mapping; }
     41  static short   SetDefaultVectorType(short vt=ColumnVector);
     42  static inline short GetDefaultVectorType() { return default_vector_type; }
    3243 
    3344  // Creation / destruction
    … …  
    5364  uint_4 MaxSizeKA() const ;
    5465
    55   // memory organisation - packing information
     66  // memory organisation
    5667  inline short  GetMemoryMapping() const
    57                 {return ( (marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping) ; }
     68                { return ( (marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping) ; }
    5869
    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
     70  inline uint_4 RowsKA() const {return marowi_; }    // Dimension des index de lignes
     71  inline uint_4 ColsKA() const {return macoli_; }    // Dimension des index de colonnes
     72  inline uint_4 VectKA() const {return veceli_; }    // Dimension des index des elts d'un vecteur
     73  void          SetMemoryMapping(short mm=AutoMemoryMapping);
    6274
     75  // Vector type   Line or Column vector
     76  inline short  GetVectorType() const
     77                { return((marowi_ == veceli_) ? ColumnVector : RowVector); }
     78  void          SetVectorType(short vt=AutoVectorType);
     79
     80  // memory organisation  - packing information
    6381  inline bool   IsPacked() const { return(moystep_ == 1); }
    6482  inline bool   IsPackedX() const { return(step_[0] == 1); }
    … …  
    7694  uint_4 MinStepKA() const ;
    7795
     96  // Offset of element ip
    7897  uint_8 Offset(uint_8 ip=0) const ;
     98  // Offset of the i'th vector along axe ka
     99  uint_8  Offset(uint_4 ka, uint_8 i) const ;
    79100  inline uint_8  Offset(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const; 
    80101
    … …  
    85106  void           Show(ostream& os, bool si=false) const;
    86107  inline void    Show() const { Show(cout); }
    87   virtual string DataType() const = 0;
     108  virtual string InfoString() const;
    88109
    89110//  Objet DVList info
    … …  
    102123  //  Organisation memoire
    103124  static short SelectMemoryMapping(short mm);
    104   void UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny);
     125  static short SelectVectorType(short vt);
     126  void UpdateMemoryMapping(short mm);
    105127  void UpdateMemoryMapping(BaseArray const & a, short mm);
    106128
    … …  
    113135  uint_4 step_[BASEARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension
    114136  uint_4 minstep_;                   // minimal step (in any axes)
    115   uint_4 moystep_;                   // mean step 0 non regular steps
     137  uint_4 moystep_;                   // mean step if == 0 --> non regular steps
    116138  uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
    117139  uint_4 marowi_, macoli_;     // For matrices, Row index and column index in dimensions
    118140  uint_4 veceli_;              // For vectors,  dimension index = marowi_/macoli_ (Row/Col vectors)
    119 
     141  bool ck_memo_vt_;            // if true, check MemoryOrg./VectorType for CompareSize
    120142  DVList* mInfo;               // Infos (variables) attachees au tableau
    121143
    122144  static char * ck_op_msg_[6];  // Operation messages for CheckDI() CheckBound()
    123145  static uint_4 max_nprt_;      // Nb maxi d'elements imprimes
     146  static uint_4 prt_lev_;        // Niveau de print 0 ou 1
    124147  static short default_memory_mapping;  // Default memory mapping
     148  static short default_vector_type;     // Default vector type Row/Column
     149  static uint_8 openmp_size_threshold;  // Size limit for parallel routine calls
    125150};
    126151
    … …  
    132157  if ( (ka < 0) || (ka >= ndim_) ) {
    133158    string txt = "BaseArray::CheckDimensionIndex/Error ";   txt += ck_op_msg_[msg];
    134     throw(ParmError(txt));
     159    throw(RangeCheckError(txt));
    135160  }
    136161  return(ka);
    … …  
    142167       (it >= size_[3]) || (iu >= size_[4]) ) {
    143168    string txt = "BaseArray::CheckArrayBound/Error ";   txt += ck_op_msg_[msg];
    144     throw(ParmError(txt));
     169    throw(RangeCheckError(txt));
    145170  }
    146171  return;
    147172}
     173
    148174
    149175

  • trunk/SophyaLib/TArray/matharr.cc

    r804 r813  
    2727    uint_8 step = a.Step(ka);
    2828    uint_8 gpas = a.Size(ka)*step;
    29     for(j=0; j<a.Size(); j += a.Size(ka))  {
    30       pe = a.DataBlock().Begin()+a.Offset(j);
     29    uint_8 naxa = a.Size()/a.Size(ka);
     30    for(j=0; j<naxa; j++)  {
     31      pe = a.DataBlock().Begin()+a.Offset(ka,j);
    3132      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
    3233    }
    … …  
    5253    uint_8 step = a.Step(ka);
    5354    uint_8 gpas = a.Size(ka)*step;
    54     for(j=0; j<a.Size(); j += a.Size(ka))  {
    55       pe = a.DataBlock().Begin()+a.Offset(j);
     55    uint_8 naxa = a.Size()/a.Size(ka);
     56    for(j=0; j<naxa; j++)  {
     57      pe = a.DataBlock().Begin()+a.Offset(ka,j);
    5658      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
    5759    }
    … …  
    102104    uint_8 step = a.Step(ka);
    103105    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    uint_8 naxa = a.Size()/a.Size(ka);
     107    for(j=0; j<naxa; j++)  {
     108      pe = a.DataBlock().Begin()+a.Offset(ka,j);
    106109      for(k=0; k<gpas; k+=step) {
    107110        valok = (double) pe[k];

  • trunk/SophyaLib/TArray/sopemtx.cc

    r804 r813  
    6363
    6464  uint_4 IMaxAbs(uint_4 first=0);
     65  void Print(ostream & os) const ;
    6566
    6667  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
    … …  
    105106template <class T>
    106107inline 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());
     108  { switch (rckind) { case TmatrixRow  : return m.Step(m.ColsKA());
     109                      case TmatrixCol  : return m.Step(m.RowsKA());
    109110                      case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
    110111    return 0; }
    … …  
    333334
    334335template <class T>
     336void TMatrixRC<T>::Print(ostream & os) const
     337{
     338  os << " TMatrixRC<T>::Print(ostream & os) " << NElts() << " Kind="
     339     << kind << " Index=" << index << " Step= " << step << endl;
     340  for(uint_4 i=0; i<NElts(); i++) {
     341    os << (*this)(i);
     342    if (kind == TmatrixCol)  os << endl;
     343    else os << ", ";
     344  }
     345  os << endl;
     346}
     347
     348template <class T>
    335349void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
    336350{
    … …  
    435449
    436450for(uint_4 l=0; l<n; l++) {
    437   if (fabs(a(l,l)) <= 1.e-50) return 0.0;
     451  if (fabs((double)a(l,l)) <= 1.e-50) return 0.0;
    438452  TMatrixRC<T>::Row(b, l) /= a(l,l);
    439453}
    … …  
    447461{
    448462TMatrix<T> a(A);
    449 TMatrix<T> b(a.NCols(),a.NRows());  b = 1.;
    450 if(fabs(GausPiv(a,b)) < 1.e-50)
     463TMatrix<T> b(a.NCols(),a.NRows());  b = IdentityMatrix(1.);
     464if(fabs((double)GausPiv(a,b)) < 1.e-50)
    451465  throw(MathExc("TMatrix Inverse() Singular OMatrix"));
    452466return b;
    … …  
    578592
    579593
    580 
     594void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m)
     595{
     596  cout << " ::::: _ZZ_TestTMatrixRC_YY_ :::: M= \n" << m << endl;
     597  TMatrixRC<r_8> l0 = TMatrixRC<r_8>::Row(m,0);
     598  cout << "TMatrixRC<r_8>::Row(m,0) = \n" ;
     599  l0.Print(cout);
     600  TMatrixRC<r_8> l1 = TMatrixRC<r_8>::Row(m,1);
     601  cout << "TMatrixRC<r_8>::Row(m,1) = \n" ;
     602  l1.Print(cout);
     603  l0.SetRow(2);
     604  cout << "TMatrixRC<r_8>::l0.SetRow(2 = \n" ;
     605  l0.Print(cout);
     606
     607  TMatrixRC<r_8> c0 = TMatrixRC<r_8>::Col(m,0);
     608  cout << "TMatrixRC<r_8>::Col(m,0) = \n" ;
     609  c0.Print(cout);
     610  TMatrixRC<r_8> c1 = TMatrixRC<r_8>::Col(m,1);
     611  cout << "TMatrixRC<r_8>::Col(m,1) = \n" ;
     612  c1.Print(cout); 
     613  c0.SetCol(2);
     614  cout << "TMatrixRC<r_8>::c0.SetCol(2) = \n" ;
     615  c0.Print(cout);
     616  TMatrixRC<r_8> c00 = TMatrixRC<r_8>::Col(m,0);
     617  TMatrixRC<r_8>::Swap(c0, c00);
     618  cout << " :::::  M Apres Swap = \n" << m << endl;
     619 
     620}
    581621
    582622///////////////////////////////////////////////////////////////

  • trunk/SophyaLib/TArray/tarray.cc

    r804 r813  
    234234//  Possibilite de // , vectorisation
    235235template <class T>
    236 TArray<T>& TArray<T>::Set(Sequence seq)
    237 {
    238   if (NbDimensions() < 1)
    239     throw RangeCheckError("TArray<T>::Set(Sequence ) - Not Allocated Array ! ");
     236TArray<T>& TArray<T>::SetSeq(Sequence seq)
     237{
     238  if (NbDimensions() < 1)
     239    throw RangeCheckError("TArray<T>::SetSeq(Sequence ) - Not Allocated Array ! ");
    240240  T * pe;
    241241  uint_8 j,k;
    … …  
    246246  }
    247247  else {    // Non regular data spacing ...
    248     uint_4 ka = MaxSizeKA();
     248    //    uint_4 ka = MaxSizeKA();
     249    uint_4 ka = 0;
    249250    uint_8 step = Step(ka);
    250251    uint_8 gpas = Size(ka);
    251     for(j=0; j<totsize_; j += Size(ka))  {
    252       pe = mNDBlock.Begin()+Offset(j);
    253       for(k=0; k<gpas; k++)  pe[k*step] = seq(j+k);
     252    uint_8 naxa = Size()/Size(ka);
     253    for(j=0; j<naxa; j++)  {
     254      pe = mNDBlock.Begin()+Offset(ka,j);
     255      for(k=0; k<gpas; k++)  pe[k*step] = seq(j*gpas+k);
    254256    }
    255257  }
    … …  
    260262
    261263template <class T>
    262 TArray<T>& TArray<T>::Set(T x)
    263 {
    264   if (NbDimensions() < 1)
    265     throw RangeCheckError("TArray<T>::Set(T )  - Not Allocated Array ! ");
     264TArray<T>& TArray<T>::SetT(T x)
     265{
     266  if (NbDimensions() < 1)
     267    throw RangeCheckError("TArray<T>::SetT(T )  - Not Allocated Array ! ");
    266268  T * pe;
    267269  uint_8 j,k;
    … …  
    276278    uint_8 step = Step(ka);
    277279    uint_8 gpas = Size(ka)*step;
    278     for(j=0; j<totsize_; j += Size(ka))  {
    279       pe = mNDBlock.Begin()+Offset(j);
     280    uint_8 naxa = Size()/Size(ka);
     281    for(j=0; j<naxa; j++)  {
     282      pe = mNDBlock.Begin()+Offset(ka,j);
    280283      for(k=0; k<gpas; k+=step)  pe[k] = x;
    281284    }
    … …  
    301304    uint_8 step = Step(ka);
    302305    uint_8 gpas = Size(ka)*step;
    303     for(j=0; j<totsize_; j += Size(ka))  {
    304       pe = mNDBlock.Begin()+Offset(j);
     306    uint_8 naxa = Size()/Size(ka);
     307    for(j=0; j<naxa; j++)  {
     308      pe = mNDBlock.Begin()+Offset(ka,j);
    305309      for(k=0; k<gpas; k+=step)  pe[k] += x;
    306310    }
    … …  
    326330    uint_8 step = Step(ka);
    327331    uint_8 gpas = Size(ka)*step;
    328     for(j=0; j<totsize_; j += Size(ka))  {
    329       pe = mNDBlock.Begin()+Offset(j);
     332    uint_8 naxa = Size()/Size(ka);
     333    for(j=0; j<naxa; j++)  {
     334      pe = mNDBlock.Begin()+Offset(ka,j);
    330335      for(k=0; k<gpas; k+=step)  pe[k] -= x;
    331336    }
    … …  
    351356    uint_8 step = Step(ka);
    352357    uint_8 gpas = Size(ka)*step;
    353     for(j=0; j<totsize_; j += Size(ka))  {
    354       pe = mNDBlock.Begin()+Offset(j);
     358    uint_8 naxa = Size()/Size(ka);
     359    for(j=0; j<naxa; j++)  {
     360      pe = mNDBlock.Begin()+Offset(ka,j);
    355361      for(k=0; k<gpas; k+=step)  pe[k] *= x;
    356362    }
    … …  
    376382    uint_8 step = Step(ka);
    377383    uint_8 gpas = Size(ka)*step;
    378     for(j=0; j<totsize_; j += Size(ka))  {
    379       pe = mNDBlock.Begin()+Offset(j);
     384    uint_8 naxa = Size()/Size(ka);
     385    for(j=0; j<naxa; j++)  {
     386      pe = mNDBlock.Begin()+Offset(ka,j);
    380387      for(k=0; k<gpas; k+=step)  pe[k] /= x;
    381388    }
    … …  
    402409    uint_8 step = Step(ka);
    403410    uint_8 gpas = Size(ka)*step;
    404     for(j=0; j<totsize_; j += Size(ka))  {
    405       pe = mNDBlock.Begin()+Offset(j);
     411    uint_8 naxa = Size()/Size(ka);
     412    for(j=0; j<naxa; j++)  {
     413      pe = mNDBlock.Begin()+Offset(ka,j);
    406414      for(k=0; k<gpas; k+=step)  pe[k] = x-pe[k];
    407415    }
    … …  
    427435    uint_8 step = Step(ka);
    428436    uint_8 gpas = Size(ka)*step;
    429     for(j=0; j<totsize_; j += Size(ka))  {
    430       pe = mNDBlock.Begin()+Offset(j);
     437    uint_8 naxa = Size()/Size(ka);
     438    for(j=0; j<naxa; j++)  {
     439      pe = mNDBlock.Begin()+Offset(ka,j);
    431440      for(k=0; k<gpas; k+=step)  pe[k] = x/pe[k];
    432441    }
    … …  
    442451  if (NbDimensions() < 1)
    443452    throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
    444   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&)")) ;
     453  if (!CompareSizes(a))
     454    throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&) SizeMismatch")) ;
    445455
    446456  T * pe;
    … …  
    460470    uint_8 stepa = a.Step(ax);
    461471    uint_8 gpas = Size(ax)*step;
    462     for(j=0; j<totsize_; j += Size(ax))  {
    463       pe = mNDBlock.Begin()+Offset(j);
    464       pea = a.DataBlock().Begin()+a.Offset(j);
     472    uint_8 naxa = Size()/Size(ax);
     473    for(j=0; j<naxa; j++)  {
     474      pe = mNDBlock.Begin()+Offset(ax,j);
     475      pea = a.DataBlock().Begin()+a.Offset(ax,j);
    465476      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka];
    466477    }
    … …  
    474485  if (NbDimensions() < 1)
    475486    throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
    476   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&)")) ;
     487  if (!CompareSizes(a))
     488    throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&) SizeMismatch")) ;
    477489
    478490  T * pe;
    … …  
    492504    uint_8 stepa = a.Step(ax);
    493505    uint_8 gpas = Size(ax)*step;
    494     for(j=0; j<totsize_; j += Size(ax))  {
    495       pe = mNDBlock.Begin()+Offset(j);
    496       pea = a.DataBlock().Begin()+a.Offset(j);
     506    uint_8 naxa = Size()/Size(ax);
     507    for(j=0; j<naxa; j++)  {
     508      pe = mNDBlock.Begin()+Offset(ax,j);
     509      pea = a.DataBlock().Begin()+a.Offset(ax,j);
    497510      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka];
    498511    }
    … …  
    507520  if (NbDimensions() < 1)
    508521    throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
    509   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&)")) ;
     522  if (!CompareSizes(a))
     523    throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&) SizeMismatch")) ;
    510524
    511525  T * pe;
    … …  
    525539    uint_8 stepa = a.Step(ax);
    526540    uint_8 gpas = Size(ax)*step;
    527     for(j=0; j<totsize_; j += Size(ax))  {
    528       pe = mNDBlock.Begin()+Offset(j);
    529       pea = a.DataBlock().Begin()+a.Offset(j);
     541    uint_8 naxa = Size()/Size(ax);
     542    for(j=0; j<naxa; j++)  {
     543      pe = mNDBlock.Begin()+Offset(ax,j);
     544      pea = a.DataBlock().Begin()+a.Offset(ax,j);
    530545      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka];
    531546    }
    … …  
    540555  if (NbDimensions() < 1)
    541556    throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
    542   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
     557  if (!CompareSizes(a))
     558    throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&) SizeMismatch")) ;
    543559
    544560  T * pe;
    … …  
    558574    uint_8 stepa = a.Step(ax);
    559575    uint_8 gpas = Size(ax)*step;
    560     for(j=0; j<totsize_; j += Size(ax))  {
    561       pe = mNDBlock.Begin()+Offset(j);
    562       pea = a.DataBlock().Begin()+a.Offset(j);
     576    uint_8 naxa = Size()/Size(ax);
     577    for(j=0; j<naxa; j++)  {
     578      pe = mNDBlock.Begin()+Offset(ax,j);
     579      pea = a.DataBlock().Begin()+a.Offset(ax,j);
    563580      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka];
    564581    }
    … …  
    572589  if (NbDimensions() < 1)
    573590    throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
    574   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
     591  if (!CompareSizes(a))
     592    throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&) SizeMismatch")) ;
    575593
    576594  T * pe;
    … …  
    590608    uint_8 stepa = a.Step(ax);
    591609    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);
     610    uint_8 naxa = Size()/Size(ax);
     611    for(j=0; j<naxa; j++)  {
     612      pe = mNDBlock.Begin()+Offset(ax,j);
     613      pea = a.DataBlock().Begin()+a.Offset(ax,j);
    595614      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka];
    596615    }
    … …  
    619638    uint_8 step = Step(ka);
    620639    uint_8 gpas = Size(ka)*step;
    621     for(j=0; j<totsize_; j += Size(ka))  {
    622       pe = mNDBlock.Begin()+Offset(j);
     640    uint_8 naxa = Size()/Size(ka);
     641    for(j=0; j<naxa; j++)  {
     642      pe = mNDBlock.Begin()+Offset(ka,j);
    623643      for(k=0; k<gpas; k+=step)  ret += pe[k] ;
    624644    }
    … …  
    645665    uint_8 step = Step(ka);
    646666    uint_8 gpas = Size(ka)*step;
    647     for(j=0; j<totsize_; j += Size(ka))  {
    648       pe = mNDBlock.Begin()+Offset(j);
     667    uint_8 naxa = Size()/Size(ka);
     668    for(j=0; j<naxa; j++)  {
     669      pe = mNDBlock.Begin()+Offset(ka,j);
    649670      for(k=0; k<gpas; k+=step)  ret *= pe[k] ;
    650671    }
    … …  
    660681
    661682template <class T>
    662 string TArray<T>::DataType() const
    663 {
    664   string rs = typeid(T).name();
     683string TArray<T>::InfoString() const
     684{
     685  string rs = "TArray<" ;
     686  rs += typeid(T).name();
     687  rs += "> ";
    665688  return(rs);
    666689}
    … …  
    693716    }
    694717  }
    695   os << "\n" << endl;
     718  os << endl;
    696719}
    697720

  • trunk/SophyaLib/TArray/tarray.h

    r804 r813  
    8888  inline T toScalar();
    8989// Met les elements a une suite de valeurs
    90   virtual TArray<T>&  Set(Sequence seq);
    91   inline  TArray<T>&  operator = (Sequence seq)    { return Set(seq); }
     90  virtual TArray<T>&  SetSeq(Sequence seq);
     91  inline  TArray<T>&  operator = (Sequence seq)    { return SetSeq(seq); }
    9292// A = x (tous les elements a x)
    93   virtual TArray<T>&  Set(T x);
    94   inline  TArray<T>&  operator = (T x)             { return Set(x); }
     93  virtual TArray<T>&  SetT(T x);
     94  inline  TArray<T>&  operator = (T x)             { return SetT(x); }
    9595// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
    9696  virtual TArray<T>&  Add(T x);
    … …  
    121121
    122122// Impression, I/O, ...
    123   virtual string DataType() const;   // definition of abstract method inherited from BaseArray
    124   virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
     123  virtual string InfoString() const;   
     124  virtual void   Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
    125125
    126126//  Pour la gestion de persistance

  • trunk/SophyaLib/TArray/tarrinit.cc

    r804 r813  
    2121  PPRegister(FIO_TArray<uint_2>);
    2222  DObjRegister(FIO_TArray<uint_2>, TArray<uint_2>);
     23  DObjRegister(FIO_TArray<uint_2>, TMatrix<uint_2>);
     24  DObjRegister(FIO_TArray<uint_2>, TVector<uint_2>);
     25
    2326  //  PPRegister(FIO_TArray<int_2>);
    2427  //  DObjRegister(FIO_TArray<int_2>, TArray<int_2>);
     28
    2529  PPRegister(FIO_TArray<int_4>);
    2630  DObjRegister(FIO_TArray<int_4>, TArray<int_4>);
     31  DObjRegister(FIO_TArray<int_4>, TMatrix<int_4>);
     32  DObjRegister(FIO_TArray<int_4>, TVector<int_4>);
     33
    2734  PPRegister(FIO_TArray<int_8>);
    2835  DObjRegister(FIO_TArray<int_8>, TArray<int_8>);
     36  DObjRegister(FIO_TArray<int_8>, TMatrix<int_8>);
     37  DObjRegister(FIO_TArray<int_8>, TVector<int_8>);
     38
    2939  //  PPRegister(FIO_TArray<uint_4>);
    3040  // DObjRegister(FIO_TArray<uint_4>, TArray<uint_4>);
    … …  
    3343  PPRegister(FIO_TArray<r_4>);
    3444  DObjRegister(FIO_TArray<r_4>, TArray<r_4>);
     45  DObjRegister(FIO_TArray<r_4>, TMatrix<r_4>);
     46  DObjRegister(FIO_TArray<r_4>, TVector<r_4>);
     47
    3548  PPRegister(FIO_TArray<r_8>);
    3649  DObjRegister(FIO_TArray<r_8>, TArray<r_8>);
     50  DObjRegister(FIO_TArray<r_8>, TMatrix<r_8>);
     51  DObjRegister(FIO_TArray<r_8>, TVector<r_8>);
     52
    3753  PPRegister(FIO_TArray< complex<r_4> >);
    3854  DObjRegister(FIO_TArray< complex<r_4> >, TArray< complex<r_4> >);
     55  DObjRegister(FIO_TArray< complex<r_4> >, TMatrix< complex<r_4> >);
     56  DObjRegister(FIO_TArray< complex<r_4> >, TVector< complex<r_4> >);
     57
    3958  PPRegister(FIO_TArray< complex<r_8> >);
    4059  DObjRegister(FIO_TArray< complex<r_8> >, TArray< complex<r_8> >);
     60  DObjRegister(FIO_TArray< complex<r_8> >, TMatrix< complex<r_8> >);
     61  DObjRegister(FIO_TArray< complex<r_8> >, TVector< complex<r_8> >);
    4162
    4263}

  • trunk/SophyaLib/TArray/tmatrix.cc

    r804 r813  
    1 // $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
     1// $Id: tmatrix.cc,v 1.4 2000-04-05 15:44:15 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    … …  
    1717  : TArray<T>()
    1818{
     19  ck_memo_vt_ = true;
    1920}
    2021
    … …  
    4849: TArray<T>(a)
    4950{
    50   if (a.NbDimensions() != 2)
    51     throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
     51  if (a.NbDimensions() > 2)
     52    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions()>2 ");
     53  if (a.NbDimensions() == 1) {
     54    size_[1] = 1;
     55    step_[1] = size_[0]*step_[0];
     56    ndim_ = 2;
     57  }
     58  UpdateMemoryMapping(a, SameMemoryMapping);
    5259}
    5360
    … …  
    5663: TArray<T>(a, share)
    5764{
    58   if (a.NbDimensions() != 2)
    59     throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
     65  if (a.NbDimensions() > 2)
     66    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions()>2");
     67  if (a.NbDimensions() == 1) {
     68    size_[1] = 1;
     69    step_[1] = size_[0]*step_[0];
     70    ndim_ = 2;
     71  }
    6072  UpdateMemoryMapping(a, mm);
    6173}
    … …  
    7183TArray<T>& TMatrix<T>::Set(const TArray<T>& a)
    7284{
    73   if (a.NbDimensions() != 2)
    74     throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() != 2 ");
    75   return(TArray<T>::Set(a));
     85  if (a.NbDimensions() > 2)
     86    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() > 2");
     87  TArray<T>::Set(a);
     88  if (a.NbDimensions() == 1) {
     89    size_[1] = 1;
     90    step_[1] = size_[0]*step_[0];
     91    ndim_ = 2;
     92  }
     93  UpdateMemoryMapping(a, SameMemoryMapping);
     94  return(*this);
    7695}
    7796
    … …  
    83102  uint_4 size[BASEARRAY_MAXNDIMS];
    84103  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
    85   size[0] = r;  size[1] = c;
    86104  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
    87   UpdateMemoryMapping(mm, size[0], size[1]);
     105  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) )
     106    mm = GetDefaultMemoryMapping();
     107  if (mm == CMemoryMapping) {
     108    size[0] = c;  size[1] = r;
     109  }
     110  else {
     111    size[0] = r;  size[1] = c;
     112  }
    88113  TArray<T>::ReSize(2, size, 1);
    89   UpdateMemoryMapping(mm, size[0], size[1]);
     114  UpdateMemoryMapping(mm);
    90115}
    91116
    … …  
    97122  uint_4 size[BASEARRAY_MAXNDIMS];
    98123  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
    99   UpdateMemoryMapping(mm, size[0], size[1]);
     124  if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     125  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) )
     126    mm = GetDefaultMemoryMapping();
     127  if (mm == CMemoryMapping) {
     128    size[0] = c;  size[1] = r;
     129  }
     130  else {
     131    size[0] = r;  size[1] = c;
     132  }
    100133  TArray<T>::Realloc(2, size, 1, force);
    101   UpdateMemoryMapping(mm, size[0], size[1]);
     134  UpdateMemoryMapping(mm);
    102135}
    103136
    104137// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
    105138template <class T>
    106 TMatrix<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   }
     139TMatrix<T> TMatrix<T>::SubMatrix(Range rline, Range rcol) const
     140{
     141  short mm = GetMemoryMapping();
     142  Range rx, ry;
     143  if (mm == CMemoryMapping)  { rx = rcol;  ry = rline; }
     144  else { ry = rcol;  rx = rline; }
     145  TMatrix sm(SubArray(rx, ry, Range(0), Range(0), Range(0)),true, mm);
     146  sm.UpdateMemoryMapping(mm);
     147  sm.SetTemp(true);
     148  return(sm);
    122149}
    123150
    … …  
    127154TMatrix<T>& TMatrix<T>::Transpose()
    128155{
     156  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
    129157  uint_4 rci = macoli_;
    130158  macoli_ = marowi_;
    131159  marowi_ = rci;
     160  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
    132161  return(*this);
    133162}
    … …  
    147176
    148177// Rearrangement memoire
     178// Si identique, renvoie une matrice qui partage les donnees
    149179template <class T>
    150180TMatrix<T> TMatrix<T>::Rearrange(short mm)
    151181{
    152   if (mm == SameMemoryMapping) mm = GetMemoryMapping(); 
     182  if ( mm == SameMemoryMapping)  mm = GetMemoryMapping();
     183  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) )
     184    mm = GetDefaultMemoryMapping();
     185 
     186  if  (mm == GetMemoryMapping())
     187    return (TMatrix<T>(*this, true));
     188 
    153189  TMatrix<T> tm(NRows(), NCols(), mm);
    154190  for(uint_4 i=0; i<NRows(); i++)
    … …  
    180216//**** Impression
    181217template <class T>
     218string TMatrix<T>::InfoString() const
     219{
     220  string rs = "TMatrix<";
     221  rs += typeid(T).name();
     222  char buff[64];
     223  sprintf(buff, ">(NRows=%ld, NCols=%ld)", (long)NRows(), (long)NCols());
     224  rs += buff;
     225  return(rs); 
     226}
     227
     228template <class T>
    182229void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
    183230{
    184   os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols()
    185      << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
    186   os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
    187231  if (maxprt < 0)  maxprt = max_nprt_;
    188232  int_4 npr = 0;
    … …  
    200244    os << endl;
    201245  }
    202   os << "\n" << endl;
     246  os << endl;
    203247}
    204248
    … …  
    235279}
    236280
    237 
    238281///////////////////////////////////////////////////////////////
    239282#ifdef __CXX_PRAGMA_TEMPLATES__

  • trunk/SophyaLib/TArray/tmatrix.h

    r804 r813  
    3434  void Realloc(uint_4 r,uint_4 c, short mm=SameMemoryMapping, bool force=false);
    3535
    36   // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
    37   TMatrix<T> operator () (Range rline, Range rcol) const ;
     36  // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer ces methode const ?
     37  TMatrix<T> SubMatrix(Range rline, Range rcol) const ;
     38  inline TMatrix<T> operator () (Range rline, Range rcol) const
     39                    { return SubMatrix(rline, rcol); }
     40  // Lignes et colonnes de la matrice
     41  inline TMatrix<T> Row(uint_4 ir) const
     42                    { return SubMatrix(Range(ir,ir), Range(0,NCols()-1)); }
     43  inline TMatrix<T> Column(uint_4 ic) const
     44                    { return SubMatrix(Range(0,NRows()-1), Range(ic,ic)); }
    3845
    3946  // Inline element acces methods
    … …  
    5360  inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }
    5461
     62  inline  TMatrix<T>&  operator = (Sequence seq)    { SetSeq(seq); return(*this); }
     63
    5564  // Operations diverses  avec une constante
    56   inline  TMatrix<T>&  operator = (T x)             { Set(x); return(*this); }
     65  inline  TMatrix<T>&  operator = (T x)             { SetT(x); return(*this); }
    5766  inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
    5867  inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
    … …  
    6170
    6271  //  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); }
     72  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { AddElt(a); return(*this); }
     73  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { SubElt(a); return(*this); }
     74  // Produit matriciel Multiply : C = (*this)*B
     75  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
     76  inline  TMatrix<T>&  operator *= (const TMatrix<T>& b)
     77          { this->Set(Multiply(b));  return(*this); }
    6578
    66   // Input-Output
     79  // I/O print, ...
     80  virtual string InfoString() const;
    6781  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
    68 
    69   // Produit matriciel
    70   TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
    71   //  inline TMatrix<T>& operator *= (const TMatrix<T>& b)
    7282
    7383protected:
    … …  
    98108
    99109
     110// Surcharge d'operateurs C = A (+,-) B
     111// $CHECK$ Reza 3/4/2000 Pas necessaire  de redefinir les operateurs
     112// Defini au niveau de TArray<T> - Pour ameliorer l'efficacite
     113// Doit-on le faire aussi pour les constantes ? - Fin de $CHECK$ Reza 3/4/2000
     114
     115template <class T>
     116inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
     117    {TMatrix<T> result(a); result.SetTemp(true); result.AddElt(b); return result;}
     118
     119template <class T>
     120inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
     121    {TMatrix<T> result(a); result.SetTemp(true); result.SubElt(b); return result;}
     122
    100123// Surcharge d'operateurs C = A * B
    101124
    102125template <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;}
     126{ TMatrix<T> result(a); result.SetTemp(true); return(result.Multiply(b)); }
    104127
    105128typedef TMatrix<r_8> Matrix;

  • trunk/SophyaLib/TArray/tvector.cc

    r804 r813  
    1 // $Id: tvector.cc,v 1.2 2000-04-03 17:35:59 ansari Exp $
     1// $Id: tvector.cc,v 1.3 2000-04-05 15:44:17 ansari Exp $
    22//                         C.Magneville          04/99
    33#include "machdefs.h"
    … …  
    1818TVector<T>::TVector(uint_4 n, short lcv, short mm)
    1919// Constructeur
    20   : TMatrix<T>((lcv == ColumnVector) ? n : 1,(lcv == ColumnVector) ? 1 : 0, mm)
     20  : TMatrix<T>(1,1,mm)
    2121{
    22   if (lcv == ColumnVector) veceli_ = marowi_;
    23   else veceli_ = macoli_;
     22  lcv = SelectVectorType(lcv);
     23  ReSize(n,lcv);
    2424}
    2525
    … …  
    4343: TMatrix<T>(a)
    4444{
     45  if ( (size_[0] != 1) && (size_[1] != 1) )
     46    throw SzMismatchError("TVector<T>::TVector(const TArray<T>& a) NRows()!=1 && NCols()!=1 ");
    4547}
    4648
    4749
    4850template <class T>
    49 TVector<T>::TVector(const TArray<T>& a, bool share, short lcv, short mm )
     51TVector<T>::TVector(const TArray<T>& a, bool share, short mm, short lcv )
    5052: TMatrix<T>(a, share, mm)
    5153{
    52   if (lcv == SameTypeVector) veceli_ = a.VectKA();
    53   else {
    54     if (lcv == ColumnVector) veceli_ = marowi_;
    55     else veceli_ = macoli_;
     54  if ( (size_[0] != 1) && (size_[1] != 1) )
     55    throw SzMismatchError("TVector<T>::TVector(const TArray<T>& a) NRows()!=1 && NCols()!=1 ");
     56  if ( (size_[0] == 1) && (size_[1] == 1) ) {
     57    if (lcv == SameVectorType) lcv = a.GetVectorType();
     58    if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
     59    veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
    5660  }
    5761}
    … …  
    7074    throw(SzMismatchError("TVector::ReSize() n = 0 "));
    7175  uint_4 r,c;
    72   if (lcv == SameTypeVector)  lcv = GetVectorType();
     76  if (lcv == SameVectorType)  lcv = GetVectorType();
     77  else if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
    7378  if (lcv == ColumnVector) { r = n;  c = 1; }
    7479  else { c = n; r = 1; }
    7580  TMatrix<T>::ReSize(r,c);
     81  veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
    7682}
    7783
    … …  
    8288    throw(SzMismatchError("TVector::Realloc() n = 0 "));
    8389  uint_4 r,c;
    84   if (lcv == SameTypeVector)  lcv = GetVectorType();
     90  if (lcv == SameVectorType)  lcv = GetVectorType();
     91  else if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
    8592  if (lcv == ColumnVector) { r = n;  c = 1; }
    8693  else { c = n; r = 1; }
    8794  TMatrix<T>::Realloc(r,c,SameMemoryMapping,force);
     95  veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
    8896}
    8997
    9098// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
    9199template <class T>
    92 TVector<T> TVector<T>::operator () (Range relt) const
     100TVector<T> TVector<T>::SubVector(Range relt) const
    93101{
    94102  Range rr, cr;
    … …  
    109117}
    110118
     119template <class T>
     120string TVector<T>::InfoString() const
     121{
     122  string rs = "TVector<";
     123  rs += typeid(T).name();
     124  char buff[64];
     125  sprintf(buff, ">(%ld) (nr=%ld, nc=%ld)", (long)NElts(), (long)NRows(), (long)NCols());
     126  rs += buff;
     127  return(rs); 
     128
     129}
    111130
    112131///////////////////////////////////////////////////////////////

  • trunk/SophyaLib/TArray/tvector.h

    r804 r813  
    1111class TVector : public TMatrix<T> {
    1212public:
    13   enum {SameTypeVector=0, ColumnVector=1, LineVector=2};
    1413  // Creation / destruction
    1514  TVector();
    … …  
    1817  TVector(const TVector<T>& v, bool share);
    1918  TVector(const TArray<T>& a);
    20   TVector(const TArray<T>& a,  bool share, short lcv=ColumnVector, short mm=CMemoryMapping);
     19  TVector(const TArray<T>& a,  bool share, short mm=CMemoryMapping, short lcv=ColumnVector);
    2120
    2221  virtual ~TVector();
    … …  
    2524                       { Set(a);  return(*this); }
    2625
    27   // Vector type   Line or Column vector
    28   inline short GetVectorType() const
    29                { return((marowi_ == veceli_) ? ColumnVector : LineVector); }
    3026
    3127  // Gestion taille/Remplissage
    32   void ReSize(uint_4 n, short lcv=SameTypeVector );
    33   void Realloc(uint_4 n, short lcv=SameTypeVector, bool force=false);
     28  void ReSize(uint_4 n, short lcv=SameVectorType );
     29  void Realloc(uint_4 n, short lcv=SameVectorType, bool force=false);
    3430
    3531  // Sub-Vector extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
    36   TVector<T> operator () (Range relt) const ;
     32  TVector<T> SubVector(Range relt) const ;
     33  inline TVector<T> operator () (Range relt) const
     34                    { return SubVector(relt); }
    3735
    3836  // Informations pointeur/data
    … …  
    4341  inline T&       operator()(uint_4 n);
    4442
     43  // Operateur d'affectation
     44  inline  TMatrix<T>&  operator = (Sequence seq)    { SetSeq(seq); return(*this); }
     45
    4546  // Operations diverses  avec une constante
    46   inline  TVector<T>&  operator = (T x)             { Set(x); return(*this); }
     47  inline  TVector<T>&  operator = (T x)             { SetT(x); return(*this); }
    4748  inline  TVector<T>&  operator += (T x)            { Add(x); return(*this); }
    4849  inline  TVector<T>&  operator -= (T x)            { Sub(x); return(*this); }
    … …  
    5657  // Norme(^2)
    5758  T Norm2() const ;
     59
     60  virtual string InfoString() const;
     61
    5862};
    5963

  • trunk/SophyaLib/TArray/utilarr.cc

    r804 r813  
    2020}
    2121
    22 Range::Range(uint_4 start, uint_4 size, uint_4 step)
     22Range::Range(uint_4 start, uint_4 end, uint_4 size, uint_4 step)
    2323{
    2424  start_ = start;
    25   size_ = (size > 0) ? size : 1;
    2625  step_ = (step > 0) ? step : 1;
     26  if (end > start) {  // Taille calcule automatiquement
     27    end_ = end;
     28    if (step_ > ((end_-start_)+1))  size_ = 1;
     29    else size_ = ((end-start)+1)/step_;
     30  }
     31  else {     // Taille fixee
     32    size_ = size;
     33    end_ = start_+size_*step_;
     34  }
    2735}
    2836

  • trunk/SophyaLib/TArray/utilarr.h

    r804 r813  
    2929class Range {
    3030public:
    31   explicit Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
     31  explicit Range(uint_4 start=0, uint_4 end=0, uint_4 size=1, uint_4 step=1);
    3232  inline uint_4 & Start()  {  return start_; }
    33   inline uint_4 & Size()  {  return size_; }
    34   inline uint_4 & Step()  {  return step_; }
     33  inline uint_4 & End()    {  return end_; }
     34  inline uint_4 & Size()   {  return size_; }
     35  inline uint_4 & Step()   {  return step_; }
    3536protected:
    36   uint_4 start_, size_, step_;
     37  uint_4 start_, end_, size_, step_ ;
    3738};
    3839
Note: See TracChangeset for help on using the changeset viewer.