source: Sophya/trunk/SophyaLib/TArray/basarr.cc@ 804

//      Base class for numerical arrays 
//                     R. Ansari, C.Magneville   03/2000

#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include "pexceptions.h"
#include "basarr.h"



// Variables statiques globales
char * BaseArray::ck_op_msg_[6] = {"???", "Size(int )", "IsPacked(int )", 
                                 "Stride(int )", "ElemCheckBound()", "operator()" };
uint_4 BaseArray::max_nprt_ = 50;
short  BaseArray::default_memory_mapping = CMemoryMapping;

// Methodes statiques globales
void BaseArray::SetMaxPrint(uint_4 nprt)
{
  max_nprt_ = nprt;
}


uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset)
{
  uint_8 rs = step;
  for(int k=0; k<ndim; k++) rs *= siz[k];
  return(rs+offset);
}

short BaseArray::SetDefaultMemoryMapping(short mm)
{
  default_memory_mapping = ( (mm == CMemoryMapping) ? CMemoryMapping : FortranMemoryMapping) ; 
  return default_memory_mapping;
}


short BaseArray::SelectMemoryMapping(short mm)
// si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
// sinon (mm == FortranMemoryMapping)  elements d'une colonne suivant X 
{
  if ( (mm == CMemoryMapping) || (mm == FortranMemoryMapping) )  return (mm) ; 
  else return (default_memory_mapping);
}

void BaseArray::UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny)
{
  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
  if (mm == CMemoryMapping) {
    marowi_ = veceli_ = 1;  macoli_ = 0;
    uint_4 nr = nx;
    nx = ny;   ny = nr;
  }
  else {
    marowi_ = veceli_ = 0;  macoli_ = 1;
  }  
}

void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
{
  if (mm == SameMemoryMapping) 
    mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
  if (mm == CMemoryMapping) {
    marowi_ = veceli_ = 1;  macoli_ = 0;
  }
  else {
    marowi_ = veceli_ = 0;  macoli_ = 1;
  }  
}


// -------------------------------------------------------
//                Methodes de la classe
// -------------------------------------------------------

// Les constructeurs 
BaseArray::BaseArray()
  : mInfo(NULL)
// Default constructor
{
  ndim_ = 0;
  for(int k=0; k<BASEARRAY_MAXNDIMS; k++) step_[k] = size_[k] = 0;
  totsize_ = 0;
  minstep_ = 0;
  moystep_ = 0;
  offset_ = 0;
  // Default for matrices : Lines = Y , Columns = X, Default Vectore= ColumnVector
  macoli_ = 0;
  veceli_ = marowi_ = 1;
}

// Destructeur 
BaseArray::~BaseArray()
{
}


bool BaseArray::CompareSizes(const BaseArray& a)
{
  if (ndim_ != a.ndim_)  return(false);
  for(int k=0; k<ndim_; k++) 
    if (size_[k] != a.size_[k])  return(false);
  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
       (veceli_ != a.veceli_) )  return(false);
  return(true);
}

void BaseArray::CompactAllDim()
{
  if (ndim_ < 2)  return;
  uint_4 ndim = 0;
  uint_4 size[BASEARRAY_MAXNDIMS];
  uint_4 step[BASEARRAY_MAXNDIMS];
  for(int k=0; k<ndim_; k++) {
    if (size_[k] < 2)  continue;
    size[ndim] = size_[k];
    step[ndim] = step_[k];
    ndim++;
  }
  if (ndim == 0)  {
    size[0] = size_[0];
    step[0] = step_[0];
    ndim = 1;
  }
  string exmsg = "BaseArray::CompactAllDim() ";
  if (!UpdateSizes(ndim, size, step, offset_, exmsg))  throw( ParmError(exmsg) );
  return;
}

void BaseArray::CompactTrailingDim()
{
  if (ndim_ < 2)  return;
  uint_4 ndim = 0;
  uint_4 size[BASEARRAY_MAXNDIMS];
  uint_4 step[BASEARRAY_MAXNDIMS];
  for(int k=0; k<ndim_; k++) {
    size[ndim] = size_[k];
    step[ndim] = step_[k];
    if (size_[k] > 1)  ndim=k;
  }
  if (ndim == 0)  ndim = 1;
  string exmsg = "BaseArray::CompactTrailingDim() ";
  if (!UpdateSizes(ndim, size, step, offset_, exmsg))  throw( ParmError(exmsg) );
  return;
}


uint_4 BaseArray::MinStepKA() const
{
  for(int ka=0; ka<ndim_; ka++)
    if (step_[ka] == minstep_) return(ka);
  return(0);
}

uint_4 BaseArray::MaxSizeKA() const
{
  uint_4 ka = 0;
  uint_4 mx = size_[0];
  for(int k=0; k<ndim_; k++)  
    if (size_[k] > mx) {  ka = k;  mx = size_[k];  }
  return(ka);
}


//  Acces lineaire aux elements ....  Calcul d'offset

uint_8 BaseArray::Offset(uint_8 ip) const
{
if (ip == 0)  return(offset_);
uint_4 idx[BASEARRAY_MAXNDIMS];
int k;
uint_8 rest = ip;
for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
uint_8 off = offset_;
for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
return (off);
}


// ----------------------------------------------------
//       Impression, etc ...
// ----------------------------------------------------

void BaseArray::Show(ostream& os, bool si) const
{
  os << " TArray< " << DataType() << " >  NDim= " << ndim_
     << "(" << typeid(*this).name() << " )" << endl;
  os << "TotSize=" << totsize_ << " Size(X*Y*...)= " ;
  for(int k=0; k<ndim_; k++) { 
    os << size_[k];
    if (k<ndim_-1)  os << " x ";
  }
  os << endl;
  os <<  " Step(X Y ...)= " ;
  for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
  os << " Offset= " << offset_  << "\n " << endl;
  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
 
}


DVList& BaseArray::Info()
{
if (mInfo == NULL)  mInfo = new DVList;
return(*mInfo);
}


bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
{
  if (ndim >= BASEARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }
  if (step < 1) {
    exmsg += " Step(=0) Error";  return false;
  }

  minstep_ = moystep_ = step;

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  for(k=0; k<ndim; k++) {
    size_[k] = siz[k] ;
    step_[k] = totsize_*step;
    totsize_ *= size_[k];
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }
  offset_ = offset;
  // Default for matrices : Lines = Y , Columns = X
  macoli_ = 0;
  marowi_ = veceli_ = 1;
  // Update OK 
  ndim_ = ndim;
  return true;
}

bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
{
  if (ndim >= BASEARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  uint_4 minstep = step[0];
  for(k=0; k<ndim; k++) {
    size_[k] = siz[k] ;
    step_[k] = step[k];
    totsize_ *= size_[k];
    if (step_[k] < minstep)  minstep = step_[k];
  }
  if (minstep < 1) {
    exmsg += " Step(=0) Error";  return false;
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }
  uint_8 plast = 0;
  for(k=0; k<ndim; k++)   plast += (siz[k]-1)*step[k];
  if (plast == minstep*totsize_ )  moystep_ = minstep;
  else moystep_ = 0;
  minstep_ = minstep;
  offset_ = offset;
  // Default for matrices : Lines = Y , Columns = X
  macoli_ = 0;
  marowi_ = veceli_ = 1;
  // Update OK 
  ndim_ = ndim;
  return true;
}

bool BaseArray::UpdateSizes(const BaseArray& a, string & exmsg)
{
  if (a.ndim_ >= BASEARRAY_MAXNDIMS) {
    exmsg += " NDim Error";  return false;
  }

  // Flagging bad updates ...
  ndim_ = 0;

  totsize_ = 1;
  int k;
  for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
    size_[k] = 1;
    step_[k] = 0;
  }
  uint_4 minstep = a.step_[0];
  for(k=0; k<a.ndim_; k++) {
    size_[k] = a.size_[k] ;
    step_[k] = a.step_[k];
    totsize_ *= size_[k];
    if (step_[k] < minstep)  minstep = step_[k];
  }
  if (minstep < 1) {
    exmsg += " Step(=0) Error";  return false;
  }
  if (totsize_ < 1) {
    exmsg += " Size Error";  return false;
  }

  minstep_ = a.minstep_;
  moystep_ = a.moystep_;
  offset_ = a.offset_;
  macoli_ = a.macoli_;
  marowi_ = a.marowi_;
  veceli_ = a.veceli_;
  // Update OK 
  ndim_ = a.ndim_;
  return true;
}


void BaseArray::UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const
{
  if ( (ndim > ndim_) || (ndim < 1) ) 
    throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) NDim Error") );
  int k;
  for(k=0; k<ndim; k++) 
    if ( (siz[k]*step[k]+pos[k]) > size_[k] )  
      throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) Size/Pos Error") );
  uint_8 offset = offset_;
  for(k=0; k<ndim_; k++) { 
    offset += pos[k]*step_[k]; 
    step[k] *= step_[k];
  }
  string exm = "BaseArray::UpdateSubArraySizes() ";
  if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
     throw( ParmError(exm) );
  return;
}