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

//      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;
uint_4 BaseArray::prt_lev_ = 0;
short  BaseArray::default_memory_mapping = CMemoryMapping;
short  BaseArray::default_vector_type = ColumnVector;
uint_8 BaseArray::openmp_size_threshold = 200000;

// ------ Methodes statiques globales --------

//  1/ Gestion des impressions 
//! Set maximum number of printed elements and print level
/*!
  \param nprt : maximum number of print
  \param lev  : print level
*/
void BaseArray::SetMaxPrint(uint_4 nprt, uint_4 lev)
{
  max_nprt_ = nprt;
  prt_lev_ = (lev < 3) ? lev : 3;
}

//! Set Size threshold for parallel routine call
/*!
  \param thr : thresold value
*/
void BaseArray::SetOpenMPSizeThreshold(uint_8 thr)
{
  openmp_size_threshold = thr;
}


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) ? FortranMemoryMapping : CMemoryMapping; 
  return default_memory_mapping;
}

short BaseArray::SetDefaultVectorType(short vt)
{
  default_vector_type = (vt != ColumnVector) ? RowVector : ColumnVector ; 
  return default_vector_type;
}

// Memory organisation 
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);
}
short BaseArray::SelectVectorType(short vt)
{
  if ((vt == ColumnVector) || (vt == RowVector))  return(vt);
  else return(default_vector_type);
}

void BaseArray::UpdateMemoryMapping(short mm)
{
  short vt = default_vector_type;
  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
  if (mm == CMemoryMapping) {
    marowi_  = 1;  macoli_ = 0;
  }
  else {
    marowi_ = 0;  macoli_ = 1;
  }

  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
    // Choix automatique Vecteur ligne ou colonne
    if ( size_[macoli_] == 1)  veceli_ = marowi_;
    else veceli_ = macoli_;
  }
  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
}

void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
{
  short vt = default_vector_type;
  if (mm == SameMemoryMapping) {
    mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
    vt = (a.marowi_ == a.veceli_) ? ColumnVector : RowVector;
  }
  else if (mm == AutoMemoryMapping)   mm = default_memory_mapping;

  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
  if (mm == CMemoryMapping) {
    marowi_  = 1;  macoli_ = 0;
  }
  else {
    marowi_ = 0;  macoli_ = 1;
  }
  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
    // Choix automatique Vecteur ligne ou colonne
    if ( size_[macoli_] == 1)  veceli_ = marowi_;
    else veceli_ = marowi_;
  }
  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
}

void BaseArray::SetMemoryMapping(short mm)
{
  if (mm == SameMemoryMapping) return;
  if (mm == AutoMemoryMapping)  mm = default_memory_mapping;
  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) return;
  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
  if (mm == CMemoryMapping) {
    marowi_ =  1;  macoli_ = 0;
  }
  else {
    marowi_ =  0;  macoli_ = 1;
  }
  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
    // Choix automatique Vecteur ligne ou colonne
    if ( size_[macoli_] == 1)  veceli_ = marowi_;
    else veceli_ = macoli_;
  }
  else   veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
}

void BaseArray::SetVectorType(short vt)
{
  if (vt == SameVectorType) return;
  if (vt == AutoVectorType) vt =  default_vector_type;
  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
    // Choix automatique Vecteur ligne ou colonne
    if ( size_[macoli_] == 1)  veceli_ = marowi_;
    else veceli_ = macoli_;
  }
  else  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
}

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

//! Default constructor
BaseArray::BaseArray()
  : mInfo(NULL)
{
  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 : Memory organisation and Vector type 
  if (default_memory_mapping == CMemoryMapping) {
    marowi_ =  1;  macoli_ = 0;
  }
  else {
    marowi_ =  0;  macoli_ = 1;
  }
  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
}

//! Destructor 
BaseArray::~BaseArray()
{
}


//! Returns true if \b ndim and \b sizes are equal
/*!
  \param a : array to be compared
  \return true if ndim and sizes are equal, false if not
*/
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);
  //  $CHECK$   Reza doit-on verifier ca 
  if (ck_memo_vt_ && a.ck_memo_vt_)
    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
// --------------------------------------------------
// Position de l'element 0 du vecteur i selon l'axe ka
// --------------------------------------------------
uint_8 BaseArray::Offset(uint_4 ka, uint_8 i) const
{

  if ( (ndim_ < 1) || (i == 0) )  return(offset_);
  //#ifdef SO_BOUNDCHECKING
  if (ka >= ndim_) 
    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Axe KA Error");
  if ( i*size_[ka] >= totsize_ )  
    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Index Error");
  //#endif
  uint_4 idx[BASEARRAY_MAXNDIMS];
  int k;
  uint_8 rest = i;
  idx[ka] = 0;
  for(k=0; k<ndim_; k++) {
    if (k == ka) continue;
    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);
}

uint_8 BaseArray::Offset(uint_8 ip) const
{
  if ( (ndim_ < 1) || (ip == 0) )  return(offset_);
  //#ifdef SO_BOUNDCHECKING
  if (ip >= totsize_)
    throw RangeCheckError("BaseArray::Offset(uint_8 ip) Out of range index ip");
  //#endif

  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];
  }
  //#ifdef SO_BOUNDCHECKING
  if (rest != 0) 
    throw PError("BaseArray::Offset(uint_8 ip) GUG !!! rest != 0");
  //#endif
//   if (rest != 0) cerr << " BUG ---- BaseArray::Offset( " << ip << " )" << rest << endl;
//   cerr << " DBG-Offset( " << ip << ")" ;
//   for(k=0; k<ndim_; k++) cerr << idx[k] << "," ;
//   cerr << " ZZZZ " << endl;
  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
{
  if (ndim_ < 1) {
    os << "\n--- " << BaseArray::InfoString() << " Unallocated Array ! " << endl;
    return;
  } 
  os << "\n--- " << InfoString() ; 
  os << " ND=" << ndim_ << " SizeX*Y*...= " ;
  for(int k=0; k<ndim_; k++) { 
    os << size_[k];
    if (k<ndim_-1)  os << "x";
  }
  os << " ---" << endl;
  if (prt_lev_ > 0) {
    os <<  " TotSize= " << totsize_ << " Step(X Y ...)="  ;
    for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
    os << " Offset= " << offset_  << endl;
  }
  if (prt_lev_ > 1) {
    os << " MemoryMapping=" << GetMemoryMapping() << " VecType= " << GetVectorType()
       << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() 
       << " VectKA=" << VectKA() << endl;
  }
  if (!si && (prt_lev_ < 2)) return;
  if (mInfo != NULL) os << (*mInfo) << endl;
 
}

string BaseArray::InfoString() const
{
  string rs = "BaseArray Type= ";
  rs +=  typeid(*this).name() ;
  return rs;
}

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 : Memory organisation and Vector type 
  if (default_memory_mapping == CMemoryMapping) {
    marowi_ =  1;  macoli_ = 0;
  }
  else {
    marowi_ =  0;  macoli_ = 1;
  }
  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
  // 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 : Memory organisation and Vector type 
  if (default_memory_mapping == CMemoryMapping) {
    marowi_ =  1;  macoli_ = 0;
  }
  else {
    marowi_ =  0;  macoli_ = 1;
  }
  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
  // 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_;
  ck_memo_vt_ = a.ck_memo_vt_;
  // 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;
}