Changeset 785 in Sophya


Ignore:
Timestamp:
Mar 16, 2000, 7:36:24 PM (26 years ago)
Author:
ansari
Message:

Corrections,amelioration de TArray<T> - Reza 16/3/2000

Location:
trunk/SophyaLib/TArray
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/fioarr.cc

    r772 r785  
    7474  is.Get(dobj->step_, TARRAY_MAXNDIMS);
    7575  is.Get(dobj->minstep_);
     76  is.Get(dobj->moystep_);
    7677  is.Get(dobj->offset_);
    7778  is.Get(dobj->marowi_);
    … …  
    101102  os.Put(dobj->step_, TARRAY_MAXNDIMS);
    102103  os.Put(dobj->minstep_);
     104  os.Put(dobj->moystep_);
    103105  os.Put(dobj->offset_);
    104106  os.Put(dobj->marowi_);

  • trunk/SophyaLib/TArray/tarray.cc

    r772 r785  
    99
    1010
    11 // Classe utilitaires
    12 Sequence::Sequence(double start, double step)
    13 {
    14   start_ = start;
    15   step_ = step;
    16 }
    1711
    1812// Variables statiques globales
    … …  
    5347  totsize_ = 0;
    5448  minstep_ = 0;
     49  moystep_ = 0;
    5550  offset_ = 0;
    5651  // Default for matrices : Lines = Y , Columns = X
    … …  
    123118
    124119template <class T>
    125 TArray<T> TArray<T>::SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos)
    126 {
    127   if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
    128   int k;
    129   for(k=0; k<ndim; k++)
    130     if ( (siz[k]+pos[k]) > size_[k] ) 
    131       throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (1)") );
    132   for(k=ndim; k<ndim_; k++)   
    133     if ( (pos[k]) >= size_[k] ) 
    134       throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (2)") );
    135   TArray<T> ra(*this);
    136   ra.ndim_ = ndim;
    137   for(k=ndim; k<TARRAY_MAXNDIMS; k++) {
    138     ra.size_[k] = 1;
    139     ra.step_[k] = 0;
    140   }
    141   ra.offset_ = offset_;
    142   for(k=0; k<ndim_; k++) ra.offset_ += pos[k]*step_[k];
    143   ra.SetTemp(true);
    144   return(ra);   
    145 }
    146 
    147 template <class T>
    148120TArray<T>& TArray<T>::operator = (const TArray<T>& a)
    149121{
    … …  
    192164}
    193165
     166template <class T>
     167TArray<T>& TArray<T>::CompactDim()
     168{
     169  if (ndim_ < 2)  return(*this);
     170  uint_4 ndim = 0;
     171  uint_4 size[TARRAY_MAXNDIMS];
     172  uint_4 step[TARRAY_MAXNDIMS];
     173  for(int k=0; k<ndim_; k++) {
     174    if (size_[k] < 2)  continue;
     175    size[ndim] = size_[k];
     176    step[ndim] = step_[k];
     177    ndim++;
     178  }
     179  if (ndim == 0)  {
     180    size[0] = size_[0];
     181    step[0] = step_[0];
     182    ndim = 1;
     183  }
     184  string exmsg = "TArray<T>::CompactDim() ";
     185  if (!UpdateSizes(ndim, size, step, offset_, exmsg))  throw( ParmError(exmsg) );
     186  return(*this);
     187}
     188
     189
     190template <class T>
     191uint_4 TArray<T>::MinStepKA() const
     192{
     193  for(int ka=0; ka<ndim_; ka++)
     194    if (step_[ka] == minstep_) return(ka);
     195  return(0);
     196}
     197
     198template <class T>
     199uint_4 TArray<T>::MaxSizeKA() const
     200{
     201  uint_4 ka = 0;
     202  uint_4 mx = size_[0];
     203  for(int k=0; k<ndim_; k++) 
     204    if (size_[k] > mx) {  ka = k;  mx = size_[k];  }
     205  return(ka);
     206}
     207
     208
     209//  Acces lineaire aux elements ....  Calcul d'offset
     210
     211template <class T>
     212uint_8 TArray<T>::Offset(uint_8 ip) const
     213{
     214if (ip == 0)  return(offset_);
     215uint_4 idx[TARRAY_MAXNDIMS];
     216int k;
     217uint_8 rest = ip;
     218for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
     219uint_8 off = offset_;
     220for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
     221return (off);
     222}
     223
     224// SubArrays
     225template <class T>
     226TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
     227{
     228  uint_4 ndim = 0;
     229  uint_4 size[TARRAY_MAXNDIMS];
     230  uint_4 step[TARRAY_MAXNDIMS];
     231  uint_4 pos[TARRAY_MAXNDIMS];
     232  size[0] = rx.Size();
     233  size[1] = ry.Size();
     234  size[2] = rz.Size();
     235  size[3] = rt.Size();
     236  size[4] = ru.Size();
     237
     238  step[0] = rx.Step();
     239  step[1] = ry.Step();
     240  step[2] = rz.Step();
     241  step[3] = rt.Step();
     242  step[4] = ru.Step();
     243
     244  pos[0] = rx.Start();
     245  pos[1] = ry.Start();
     246  pos[2] = rz.Start();
     247  pos[3] = rt.Start();
     248  pos[4] = ru.Start();
     249
     250  ndim = ndim_;
     251  TArray<T> ra;
     252  SubArray(ra, ndim, size, pos, step);
     253  ra.SetTemp(true);
     254  return(ra);
     255}
     256
    194257//   ...... Operation de calcul sur les tableaux ......
    195258//              ------- Attention --------
    … …  
    197260//  Possibilite de // , vectorisation
    198261template <class T>
    199 TArray<T>& TArray<T>::operator = (Sequence & seq)
    200 {
    201   T * pe = Data();
    202   uint_8 step = Step();
    203   for(uint_8 k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
     262TArray<T>& TArray<T>::operator = (Sequence seq)
     263{
     264  T * pe;
     265  uint_8 j,k;
     266  if (AvgStep() > 0)   {  // regularly spaced elements
     267    uint_8 step = AvgStep();
     268    pe = Data();
     269    for(k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
     270  }
     271  else {    // Non regular data spacing ...
     272    uint_4 ka = MaxSizeKA();
     273    uint_8 step = Step(ka);
     274    uint_8 gpas = Size(ka);
     275    for(j=0; j<totsize_; j += Size(ka))  {
     276      pe = mNDBlock.Begin()+Offset(j);
     277      for(k=0; k<gpas; k++)  pe[k*step] = seq(j+k);
     278    }
     279  }
    204280  return(*this);
    205281}
    … …  
    210286TArray<T>& TArray<T>::operator = (T x)
    211287{
    212   T * pe = Data();
    213   uint_8 step = Step();
    214   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = x;
     288  T * pe;
     289  uint_8 j,k;
     290  if (AvgStep() > 0)   {  // regularly spaced elements
     291    uint_8 step = AvgStep();
     292    uint_8 maxx = totsize_*step;
     293    pe = Data();
     294    for(k=0; k<maxx; k+=step )  pe[k] = x;
     295  }
     296  else {    // Non regular data spacing ...
     297    uint_4 ka = MaxSizeKA();
     298    uint_8 step = Step(ka);
     299    uint_8 gpas = Size(ka)*step;
     300    for(j=0; j<totsize_; j += Size(ka))  {
     301      pe = mNDBlock.Begin()+Offset(j);
     302      for(k=0; k<gpas; k+=step)  pe[k] = x;
     303    }
     304  }
    215305  return(*this);
    216306}
    … …  
    219309TArray<T>& TArray<T>::operator += (T x)
    220310{
    221   T * pe = Data();
    222   uint_8 step = Step();
    223   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] += x;
     311  T * pe;
     312  uint_8 j,k;
     313  if (AvgStep() > 0)   {  // regularly spaced elements
     314    uint_8 step = AvgStep();
     315    uint_8 maxx = totsize_*step;
     316    pe = Data();
     317    for(k=0; k<maxx; k+=step )  pe[k] += x;
     318  }
     319  else {    // Non regular data spacing ...
     320    uint_4 ka = MaxSizeKA();
     321    uint_8 step = Step(ka);
     322    uint_8 gpas = Size(ka)*step;
     323    for(j=0; j<totsize_; j += Size(ka))  {
     324      pe = mNDBlock.Begin()+Offset(j);
     325      for(k=0; k<gpas; k+=step)  pe[k] += x;
     326    }
     327  }
    224328  return(*this);
    225329}
    … …  
    228332TArray<T>& TArray<T>::operator -= (T x)
    229333{
    230   T * pe = Data();
    231   uint_8 step = Step();
    232   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] -= x;
     334  T * pe;
     335  uint_8 j,k;
     336  if (AvgStep() > 0)   {  // regularly spaced elements
     337    uint_8 step = AvgStep();
     338    uint_8 maxx = totsize_*step;
     339    pe = Data();
     340    for(k=0; k<maxx; k+=step )  pe[k] -= x;
     341  }
     342  else {    // Non regular data spacing ...
     343    uint_4 ka = MaxSizeKA();
     344    uint_8 step = Step(ka);
     345    uint_8 gpas = Size(ka)*step;
     346    for(j=0; j<totsize_; j += Size(ka))  {
     347      pe = mNDBlock.Begin()+Offset(j);
     348      for(k=0; k<gpas; k+=step)  pe[k] -= x;
     349    }
     350  }
    233351  return(*this);
    234352}
    … …  
    237355TArray<T>& TArray<T>::operator *= (T x)
    238356{
    239   T * pe = Data();
    240   uint_8 step = Step();
    241   for(uint_8 k=0; k<totsize_; k += step)  pe[k] *= x;
     357  T * pe;
     358  uint_8 j,k;
     359  if (AvgStep() > 0)   {  // regularly spaced elements
     360    uint_8 step = AvgStep();
     361    uint_8 maxx = totsize_*step;
     362    pe = Data();
     363    for(k=0; k<maxx; k+=step )  pe[k] *= x;
     364  }
     365  else {    // Non regular data spacing ...
     366    uint_4 ka = MaxSizeKA();
     367    uint_8 step = Step(ka);
     368    uint_8 gpas = Size(ka)*step;
     369    for(j=0; j<totsize_; j += Size(ka))  {
     370      pe = mNDBlock.Begin()+Offset(j);
     371      for(k=0; k<gpas; k+=step)  pe[k] *= x;
     372    }
     373  }
    242374  return(*this);
    243375}
    … …  
    246378TArray<T>& TArray<T>::operator /= (T x)
    247379{
    248   T * pe = Data();
    249   uint_8 step = Step();
    250   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] /= x;
     380  T * pe;
     381  uint_8 j,k;
     382  if (AvgStep() > 0)   {  // regularly spaced elements
     383    uint_8 step = AvgStep();
     384    uint_8 maxx = totsize_*step;
     385    pe = Data();
     386    for(k=0; k<maxx; k+=step )  pe[k] /= x;
     387  }
     388  else {    // Non regular data spacing ...
     389    uint_4 ka = MaxSizeKA();
     390    uint_8 step = Step(ka);
     391    uint_8 gpas = Size(ka)*step;
     392    for(j=0; j<totsize_; j += Size(ka))  {
     393      pe = mNDBlock.Begin()+Offset(j);
     394      for(k=0; k<gpas; k+=step)  pe[k] /= x;
     395    }
     396  }
    251397  return(*this);
    252398}
    … …  
    258404{
    259405  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator += (const TArray<T>&)")) ;
    260   T * pe = Data();
    261   const T * pea = a.Data();
    262   uint_8 step = Step();
    263   uint_8 stepa = a.Step();
    264   uint_8 k, ka;
    265   k = ka = 0;
    266   for(k=0; k<totsize_; k += step)  {
    267     pe[k] += pea[ka];
    268     ka += stepa;
     406
     407  T * pe;
     408  const T * pea;
     409  uint_8 j,k,ka;
     410  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     411    uint_8 step = AvgStep();
     412    uint_8 stepa = a.AvgStep();
     413    uint_8 maxx = totsize_*step;
     414    pe = Data();
     415    pea = a.Data();
     416    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] += pea[ka] ;
     417  }
     418  else {    // Non regular data spacing ...
     419    uint_4 ax = MaxSizeKA();
     420    uint_8 step = Step(ax);
     421    uint_8 stepa = a.Step(ax);
     422    uint_8 gpas = Size(ax)*step;
     423    for(j=0; j<totsize_; j += Size(ax))  {
     424      pe = mNDBlock.Begin()+Offset(j);
     425      pea = a.DataBlock().Begin()+a.Offset(j);
     426      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka];
     427    }
    269428  }
    270429  return(*this);
    … …  
    275434{
    276435  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
    277   T * pe = Data();
    278   const T * pea = a.Data();
    279   uint_8 step = Step();
    280   uint_8 stepa = a.Step();
    281   uint_8 k, ka;
    282   k = ka = 0;
    283   for(k=0; k<totsize_; k += step)  {
    284     pe[k] -= pea[ka];
    285     ka += stepa;
     436
     437  T * pe;
     438  const T * pea;
     439  uint_8 j,k,ka;
     440  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     441    uint_8 step = AvgStep();
     442    uint_8 stepa = a.AvgStep();
     443    uint_8 maxx = totsize_*step;
     444    pe = Data();
     445    pea = a.Data();
     446    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] -= pea[ka] ;
     447  }
     448  else {    // Non regular data spacing ...
     449    uint_4 ax = MaxSizeKA();
     450    uint_8 step = Step(ax);
     451    uint_8 stepa = a.Step(ax);
     452    uint_8 gpas = Size(ax)*step;
     453    for(j=0; j<totsize_; j += Size(ax))  {
     454      pe = mNDBlock.Begin()+Offset(j);
     455      pea = a.DataBlock().Begin()+a.Offset(j);
     456      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka];
     457    }
    286458  }
    287459  return(*this);
    … …  
    293465{
    294466  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
    295   T * pe = Data();
    296   const T * pea = a.Data();
    297   uint_8 step = Step();
    298   uint_8 stepa = a.Step();
    299   uint_8 k, ka;
    300   k = ka = 0;
    301   for(k=0; k<totsize_; k += step)  {
    302     pe[k] *= pea[ka];
    303     ka += stepa;
     467
     468  T * pe;
     469  const T * pea;
     470  uint_8 j,k,ka;
     471  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     472    uint_8 step = AvgStep();
     473    uint_8 stepa = a.AvgStep();
     474    uint_8 maxx = totsize_*step;
     475    pe = Data();
     476    pea = a.Data();
     477    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] *= pea[ka] ;
     478  }
     479  else {    // Non regular data spacing ...
     480    uint_4 ax = MaxSizeKA();
     481    uint_8 step = Step(ax);
     482    uint_8 stepa = a.Step(ax);
     483    uint_8 gpas = Size(ax)*step;
     484    for(j=0; j<totsize_; j += Size(ax))  {
     485      pe = mNDBlock.Begin()+Offset(j);
     486      pea = a.DataBlock().Begin()+a.Offset(j);
     487      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka];
     488    }
    304489  }
    305490  return(*this);
    … …  
    310495{
    311496  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
    312   T * pe = Data();
    313   const T * pea = a.Data();
    314   uint_8 step = Step();
    315   uint_8 stepa = a.Step();
    316   uint_8 k, ka;
    317   k = ka = 0;
    318   for(k=0; k<totsize_; k += step)  {
    319     pe[k] /= pea[ka];
    320     ka += stepa;
     497
     498  T * pe;
     499  const T * pea;
     500  uint_8 j,k,ka;
     501  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
     502    uint_8 step = AvgStep();
     503    uint_8 stepa = a.AvgStep();
     504    uint_8 maxx = totsize_*step;
     505    pe = Data();
     506    pea = a.Data();
     507    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] /= pea[ka] ;
     508  }
     509  else {    // Non regular data spacing ...
     510    uint_4 ax = MaxSizeKA();
     511    uint_8 step = Step(ax);
     512    uint_8 stepa = a.Step(ax);
     513    uint_8 gpas = Size(ax)*step;
     514    for(j=0; j<totsize_; j += Size(ax))  {
     515      pe = mNDBlock.Begin()+Offset(j);
     516      pea = a.DataBlock().Begin()+a.Offset(j);
     517      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka];
     518    }
    321519  }
    322520  return(*this);
    … …  
    330528TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
    331529{
    332   T * pe = Data();
    333   uint_8 step = Step();
    334   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((double)pe[k]));
     530  T * pe;
     531  uint_8 j,k;
     532  if (AvgStep() > 0)   {  // regularly spaced elements
     533    uint_8 step = AvgStep();
     534    uint_8 maxx = totsize_*step;
     535    pe = Data();
     536    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((double)pe[k]));
     537  }
     538  else {    // Non regular data spacing ...
     539    uint_4 ka = MaxSizeKA();
     540    uint_8 step = Step(ka);
     541    uint_8 gpas = Size(ka)*step;
     542    for(j=0; j<totsize_; j += Size(ka))  {
     543      pe = mNDBlock.Begin()+Offset(j);
     544      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
     545    }
     546  }
    335547  return(*this);
    336548}
    … …  
    339551TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
    340552{
    341   T * pe = Data();
    342   uint_8 step = Step();
    343   for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((float)pe[k]));
     553  T * pe;
     554  uint_8 j,k;
     555  if (AvgStep() > 0)   {  // regularly spaced elements
     556    uint_8 step = AvgStep();
     557    uint_8 maxx = totsize_*step;
     558    pe = Data();
     559    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((float)pe[k]));
     560  }
     561  else {    // Non regular data spacing ...
     562    uint_4 ka = MaxSizeKA();
     563    uint_8 step = Step(ka);
     564    uint_8 gpas = Size(ka)*step;
     565    for(j=0; j<totsize_; j += Size(ka))  {
     566      pe = mNDBlock.Begin()+Offset(j);
     567      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
     568    }
     569  }
    344570  return(*this);
    345571}
    … …  
    382608  os <<  " Step(X Y ...)= " ;
    383609  for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
    384   os << "\n " << endl;
     610  os << " Offset= " << offset_  << "\n " << endl;
    385611  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
    386612 
    … …  
    395621  int k0,k1,k2,k3,k4;
    396622  for(k4=0; k4<size_[4]; k4++) {
    397     if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4;
     623    if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
    398624    for(k3=0; k3<size_[3]; k3++) {
    399       if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3;
     625      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3 << endl;
    400626      for(k2=0; k2<size_[2]; k2++) {
    401         if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2;
     627        if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2 << endl;
    402628        for(k1=0; k1<size_[1]; k1++) {
    403           if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1;
     629          if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1 << endl;
    404630          for(k0=0; k0<size_[0]; k0++) {
    405             os << Elem(k0, k1, k2, k3, k4) << ", ";    npr++;
     631            if(k0 > 0) os << ", "; 
     632            os << Elem(k0, k1, k2, k3, k4);     npr++;
    406633            if (npr >= maxprt) {
    407634              if (npr < totsize_)  os << "\n     .... " << endl; return;
    … …  
    434661  }
    435662
    436   minstep_ = 1;
     663  minstep_ = moystep_ = step;
    437664
    438665  // Flagging bad updates ...
    … …  
    478705    step_[k] = 0;
    479706  }
    480   minstep_ = step[0];
     707  uint_4 minstep = step[0];
    481708  for(k=0; k<ndim; k++) {
    482709    size_[k] = siz[k] ;
    483710    step_[k] = step[k];
    484711    totsize_ *= size_[k];
    485     if (step_[k] < minstep_)  minstep_ = step_[k];
    486   }
    487   if (minstep_ < 1) {
     712    if (step_[k] < minstep)  minstep = step_[k];
     713  }
     714  if (minstep < 1) {
    488715    exmsg += " Step(=0) Error";  return false;
    489716  }
    … …  
    491718    exmsg += " Size Error";  return false;
    492719  }
     720  uint_8 plast = 0;
     721  for(k=0; k<ndim; k++)   plast += (siz[k]-1)*step[k];
     722  if (plast == minstep*totsize_ )  moystep_ = minstep;
     723  else moystep_ = 0;
     724  minstep_ = minstep;
    493725  offset_ = offset;
    494726  // Default for matrices : Lines = Y , Columns = X
    … …  
    499731  return true;
    500732}
     733
    501734template <class T>
    502735bool TArray<T>::UpdateSizes(const TArray<T>& a, string & exmsg)
    … …  
    515748    step_[k] = 0;
    516749  }
    517   minstep_ = a.step_[0];
     750  uint_4 minstep = a.step_[0];
    518751  for(k=0; k<a.ndim_; k++) {
    519752    size_[k] = a.size_[k] ;
    520753    step_[k] = a.step_[k];
    521754    totsize_ *= size_[k];
    522     if (step_[k] < minstep_)  minstep_ = step_[k];
    523   }
    524   if (minstep_ < 1) {
     755    if (step_[k] < minstep)  minstep = step_[k];
     756  }
     757  if (minstep < 1) {
    525758    exmsg += " Step(=0) Error";  return false;
    526759  }
    … …  
    529762  }
    530763
     764  minstep_ = a.minstep_;
     765  moystep_ = a.moystep_;
    531766  offset_ = a.offset_;
    532767  macoli_ = a.macoli_;
    … …  
    553788}
    554789
     790
     791template <class T>
     792void TArray<T>::SubArray(TArray<T> & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step)
     793{
     794  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
     795  int k;
     796  for(k=0; k<ndim; k++)
     797    if ( (siz[k]*step[k]+pos[k]) > size_[k] ) 
     798      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error") );
     799  (ra.mNDBlock).Share(mNDBlock);
     800  uint_8 offset = offset_;
     801  for(k=0; k<ndim_; k++) {
     802    offset += pos[k]*step_[k];
     803    step[k] *= step_[k];
     804  }
     805  string exm = "TArray<T>::SubArray() ";
     806  if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
     807     throw( ParmError(exm) );
     808  (ra.mNDBlock).Share(mNDBlock);
     809  return;
     810}
    555811
    556812///////////////////////////////////////////////////////////////

  • trunk/SophyaLib/TArray/tarray.h

    r772 r785  
    1414#include "ndatablock.h"
    1515#include "dvlist.h"
     16#include "utilarr.h"
    1617
    1718
    … …  
    2021
    2122namespace SOPHYA {
    22 
    23 /* Quelques utilitaires */
    24 typedef double (* Arr_DoubleFunctionOfX) (double x);
    25 typedef float  (* Arr_FloatFunctionOfX)  (float x);
    26 
    27 class Sequence {
    28 public:
    29   Sequence (double start=0., double step=1.);
    30   inline double & Start() { return start_; }
    31   inline double & Step() { return step_; }
    32   inline double operator () (uint_4 k) { return(start_+(double)k*step_); }
    33 protected:
    34   double start_, step_;
    35 };
    3623
    3724// Forward declaration
    … …  
    5643  virtual TArray<T>& operator = (const TArray<T>& a);
    5744
    58   // Extraction de sous-tableau
    59   virtual TArray<T> SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos);
    6045
    6146  // Gestion taille/Remplissage
    … …  
    6550  // Returns true if ndim and sizes are equal
    6651  virtual bool CompareSizes(const TArray<T>& a);
    67  
     52  // Compacts size=1 array dimensions
     53  virtual TArray<T>& CompactDim();
     54
    6855  // Array dimensions
    6956  inline uint_4 NbDimensions() const { return( ndim_ ); }
    … …  
    7562  inline uint_4 Size(int ka) const { return(size_[CheckDI(ka,1)]); }
    7663
     64  uint_4 MaxSizeKA() const ;
     65
    7766  // memory organisation - packing information
    78   inline bool   IsPacked() const { return(minstep_ == 1); }
     67  inline bool   IsPacked() const { return(moystep_ == 1); }
    7968  inline bool   IsPackedX() const { return(step_[0] == 1); }
    8069  inline bool   IsPackedY() const { return(step_[1] == 1); }
    … …  
    8271  inline bool   IsPacked(int ka) const { return(step_[CheckDI(ka,2)] == 1); }
    8372
    84   inline uint_4 Step() const  { return(minstep_); }
     73  inline uint_4 MinStep() const  { return(minstep_); }
     74  inline uint_4 AvgStep() const  { return(moystep_); }
    8575  inline uint_4 StepX() const { return(step_[0]); }
    8676  inline uint_4 StepZ() const { return(step_[1]); }
    … …  
    8878  inline uint_4 Step(int ka) const { return(step_[CheckDI(ka,3)]); }
    8979
    90   inline uint_8 Offset() const { return(offset_); }
     80  uint_4 MinStepKA() const ;
     81
     82  uint_8 Offset(uint_8 ip=0) const ;
    9183
    9284
    … …  
    9486  inline bool   IsTemp(void) const {return mNDBlock.IsTemp();}
    9587  inline void   SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
     88
     89  // SubArrays
     90  virtual TArray<T> operator () (Range rx, Range ry, Range rz=0, Range rt=0, Range ru=0);
    9691
    9792  // Acces to data
    … …  
    10095  inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0) const ;
    10196  inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0);
    102   inline T const& operator[](uint_4 ip) const ;
    103   inline T&       operator[](uint_4 ip);
     97  inline T const& operator[](uint_8 ip) const ;
     98  inline T&       operator[](uint_8 ip);
    10499
    105100  inline T const& Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
    … …  
    117112  inline operator T();
    118113// Met les elements a une suite de valeurs
    119   virtual TArray<T>&  operator = (Sequence & seq);
     114  virtual TArray<T>&  operator = (Sequence seq);
    120115// A = x (tous les elements a x)
    121116  virtual TArray<T>&  operator = (T x);
    … …  
    162157  // Share: partage les donnees de "a"
    163158  void Share(const TArray<T>& a);
     159  // Extraction de sous-tableau
     160  virtual void SubArray(TArray<T> & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step);
    164161
    165162  uint_4 ndim_;                   // nb of dimensions
    … …  
    168165  uint_4 step_[TARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension
    169166  uint_4 minstep_;                   // minimal step (in any axes)
     167  uint_4 moystep_;                   // mean step 0 non regular steps
    170168  uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
    171169  NDataBlock<T> mNDBlock;      // Le bloc des donnees
    … …  
    254252}
    255253
     254
    256255// --------------------------------------------------
    257256//        inline element acces methods
    … …  
    329328}
    330329
    331 template <class T>
    332 inline T const& TArray<T>::operator[](uint_4 ip) const
     330
     331template <class T>
     332inline T const& TArray<T>::operator[](uint_8 ip) const
    333333{
    334334#ifdef SO_BOUNDCHECKING
    335335  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
    336336#endif
    337 return *(mNDBlock.Begin()+offset_+ip*minstep_);
    338 }
    339 
    340 template <class T>
    341 inline T & TArray<T>::operator[](uint_4 ip)
     337return *(mNDBlock.Begin()+Offset(ip));
     338}
     339
     340template <class T>
     341inline T & TArray<T>::operator[](uint_8 ip)
    342342{
    343343#ifdef SO_BOUNDCHECKING
    344344  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
    345345#endif
    346 return *(mNDBlock.Begin()+offset_+ip*minstep_);
    347 }
     346return *(mNDBlock.Begin()+Offset(ip));
     347}
     348
    348349
    349350template <class T>
Note: See TracChangeset for help on using the changeset viewer.