arrays.h

Zobrazit dokumentaci tohoto souboru.

#ifndef GELIB_ARRAYS_H
#define GELIB_ARRAYS_H

#include "gelib.h"
#include "mathlib.h"
#include "librw.h"

#include <stdio.h>
#include <string.h>
#include <math.h>


namespace gelibspace {

#define PYP 3
  
struct VectoR;
struct PoinT;

class Dmtrx;

//  ********************************************************************************************
//  *** *** *** ***                         CLASS ELEM3D                         *** *** *** ***
//  ********************************************************************************************
struct Elem3D
{
  double x, y, z;
  
  Elem3D() { x=y=z = 0.0; }
  virtual ~Elem3D() { }
  
  //  ***  GET  ***
  double  operator[] (int i) const {
#ifdef DEBUG
    if (i<0 || i>2)  _errorr ("i out of the range");
#endif
    return ((i==0) ? x : ((i==1) ? y : z));
  }
  double& operator[] (int i)       {
#ifdef DEBUG
    if (i<0 || i>2)  _errorr ("i out of the range");
#endif
    return ((i==0) ? x : ((i==1) ? y : z));
  }
  
  //  ***  FIND OUT  ***
  bool is_identical_to (const Elem3D *p, double zero) const { return  isZero(zero, x - p->x) && isZero(zero, y - p->y) && isZero(zero, z - p->z); }
  bool is_identical_to (const Elem3D *p)              const { return        (     x == p->x  &&             y == p->y  &&             z == p->z); }
  
  //  ***  EDIT  ***
  Elem3D* tms  (double val)      { x *= val;   y *= val;   z *= val;   return this; }
  Elem3D* dvd  (double val)      { x /= val;   y /= val;   z /= val;   return this; }
  Elem3D* add  (const Elem3D *p) { x += p->x;  y += p->y;  z += p->z;  return this; }
  Elem3D* sub  (const Elem3D *p) { x -= p->x;  y -= p->y;  z -= p->z;  return this; }
  Elem3D* add  (const double *p) { x += p[0];  y += p[1];  z += p[2];  return this; }
  Elem3D* copy (const Elem3D *p) { x  = p->x;  y  = p->y;  z  = p->z;  return this; }
  Elem3D* zero (void) { x = y = z = 0.0;  return this; }
  Elem3D* round2abszero (double zr)  { if (fabs(x)<zr) x=0.0;  if (fabs(y)<zr) y=0.0;  if (fabs(z)<zr) z=0.0;  return this; }
  
  //  ***  SET  ***
  bool scan_x (FILE *stream) { return ( fscanf (stream, "%lf", &x) == 1); }
  bool scan_y (FILE *stream) { return ( fscanf (stream, "%lf", &y) == 1); }
  bool scan_z (FILE *stream) { return ( fscanf (stream, "%lf", &z) == 1); }
  bool scan_xyz (FILE *stream) { return ( fscanf (stream, "%lf %lf %lf", &x, &y, &z) == 3); }
  bool scan_xyz (const char *&src)  {     sscanf (src,    "%lf %lf %lf", &x, &y, &z);  return (SP_skip_word(src) && SP_skip_word(src) && SP_skip_word(src)); }
  bool scan_xyz (const double *src) {      x = src[0]; y = src[1]; z = src[2];  return true; }
  
  //  ***  GET  ***
  void copy_to (double *dest) const { dest[0] = x; dest[1] = y; dest[2] = z; }
  
  //  ***  COMPUTE  ***
  double giveScalProduct (const Elem3D *v) const { return x*v->x + y*v->y + z*v->z; }
  double give_sum (void) const { return x + y + z; }
  
  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force=false) const { fprintf (stdout, "PoinT [3]:\n% .*e % .*e % .*e\n", PYP, x, PYP, y, PYP, z); }
#endif
  
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS POINT                          *** *** *** ***
//  ********************************************************************************************

struct PoinT : public Elem3D
{
  //  ***  EDIT  ***
  PoinT* copy (const PoinT *p) { Elem3D::copy(p);  return this; }
  
  double dist_to  (const PoinT *p) const { return sqrt( this->dist2_to(p) ); }
  double dist2_to (const PoinT *p) const { return ((x-p->x)*(x-p->x) + (y-p->y)*(y-p->y) + (z-p->z)*(z-p->z)); }
  
  inline double dist_to_line (const PoinT *p1, const PoinT *p2) const;
  void bePointAtAbscissa (const PoinT *p1, const PoinT *p2, double ksi)
  {
    x = 0.5 * (p1->x + p2->x + ksi * ( p2->x - p1->x ));
    y = 0.5 * (p1->y + p2->y + ksi * ( p2->y - p1->y ));
    z = 0.5 * (p1->z + p2->z + ksi * ( p2->z - p1->z ));
  }

  bool give_ksiAtAbscissa (double zero, double norm, const PoinT *p1, const PoinT *p2, double &ksi) const
  {
    if ( dist_to_line (p1,p2) > norm*zero )  return false;

    long i;
    if (fabs(p1->x - p2->x) > fabs(p1->y - p2->y)) if (fabs(p1->x - p2->x) > fabs(p1->z - p2->z))  i = 0;  else  i = 2;
    else                                           if (fabs(p1->y - p2->y) > fabs(p1->z - p2->z))  i = 1;  else  i = 2;
    
    ksi = 2.0 * ( (*this)[i] - (*p1)[i] ) / ( (*p2)[i] - (*p1)[i] ) - 1.0 ;
    
    return ( fabs(ksi) < (1.0+zero) );
  }
  
  void beRotatedPoint (const Dmtrx *rot, const PoinT *point);
  void beRotatedBy (const double *r);
  void beRotatedBy2d (const double *r);
  
  double* d() const
  {
    double *a = new double[3];
    a[0] = x;    a[1] = y;    a[2] = z;
    return a;
  }
  
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS VectoR                         *** *** *** ***
//  ********************************************************************************************

struct VectoR : public Elem3D
{
  //  ***  EDIT  ***
  VectoR* copy (const VectoR *p) { Elem3D::copy(p);  return this; }
  
  VectoR* beP2P (const PoinT *p1, const PoinT *p2) {
    x = p2->x - p1->x;
    y = p2->y - p1->y;
    z = p2->z - p1->z;
    return this;
  }
  VectoR* beVectProduct (const VectoR *v1, const VectoR *v2) {
    x = v1->y * v2->z - v2->y * v1->z;
    y = v1->z * v2->x - v2->z * v1->x;
    z = v1->x * v2->y - v2->x * v1->y;
    return this;
  }
  VectoR* beVectProduct (const PoinT *p1, const PoinT *p2, const PoinT *p3)
  { // beVectProduct (p1, p2, p1, p3);
    x = p2->y*p3->z - p3->y*p2->z + p1->y*(p2->z-p3->z) - p1->z*(p2->y-p3->y);
    y = p2->z*p3->x - p3->z*p2->x + p1->z*(p2->x-p3->x) - p1->x*(p2->z-p3->z);
    z = p2->x*p3->y - p3->x*p2->y + p1->x*(p2->y-p3->y) - p1->y*(p2->x-p3->x);
    return this;
  }
  VectoR* beVectProduct (const PoinT *p1, const PoinT *p2, const PoinT *p3, const PoinT *p4)
  {
    VectoR v1, v2;
    v1.beP2P(p1, p2);
    v2.beP2P(p3, p4);

    this->beVectProduct (&v1, &v2);
    return this;
  }
  VectoR* normalize (void) { this->dvd( this->give_length() );  return this; }

  //  ***  COMPUTE  ***
  double give_length (void) const { return sqrt( x*x + y*y + z*z ); }
  double give_length_2d (void) const { return sqrt( x*x + y*y ); }
  double give_angle_2d (void) const { return atan(y/x); }
  VectoR* be_rot_by_angle_2d (const VectoR *p, double a)
  {
    double cs = cos(a);
    double sn = sin(a);
    
    x = p->x * cs - p->y * sn;
    y = p->x * sn + p->y * cs;
    z = p->z;
    return this;
  }
  VectoR* be_rotSS_by_angle_2d (const VectoR *p, double a)
  {
    double cs = cos(a);
    double sn = sin(a);
    
    x =   p->x * cs + p->y * sn;
    y = - p->x * sn + p->y * cs;
    z = p->z;
    return this;
  }
  
  //  *** INFO ***
  bool is_parallel_with (const VectoR *that, double relZero) const
  {
    double l1 = this->give_length();
    double l2 = that->give_length();
    if (l1/l2 < relZero || l2/l1 < relZero)  _warningg("Lengths of vectors are not proportional");

    double v1x = this->x/l1;  if (fabs(v1x) < relZero)  v1x = 0.0;
    double v1y = this->y/l1;  if (fabs(v1y) < relZero)  v1y = 0.0;
    double v1z = this->z/l1;  if (fabs(v1z) < relZero)  v1z = 0.0;

    double v2x = that->x/l2;  if (fabs(v2x) < relZero)  v2x = 0.0;
    double v2y = that->y/l2;  if (fabs(v2y) < relZero)  v2y = 0.0;
    double v2z = that->z/l2;  if (fabs(v2z) < relZero)  v2z = 0.0;

    int sign;
    if (v1x*v2x < 0.0 || v1y*v2y < 0.0 || v1z*v2z < 0.0)  sign = -1;
    else                                                  sign =  1;

    return (fabs(v1x-sign*v2x) < relZero  &&  fabs(v1y-sign*v2y) < relZero  &&  fabs(v1z-sign*v2z) < relZero);
  }
};

double PoinT :: dist_to_line (const PoinT *p1, const PoinT *p2) const
{
  VectoR a;
  a.beVectProduct(p1, p2, this);
  return ( sqrt((a.x*a.x + a.y*a.y + a.z*a.z) / p1->dist2_to(p2) ) );
}

//  ********************************************************************************************
//  *** *** *** ***                         CLASS MatriX                         *** *** *** ***
//  ********************************************************************************************
struct MatriX
{
  double a[9];
  
  MatriX() { }
  virtual ~MatriX() { }
  
  double& at11() { return a[0]; } 
  double& at12() { return a[1]; } 
  double& at13() { return a[2]; } 
  double& at21() { return a[3]; } 
  double& at22() { return a[4]; } 
  double& at23() { return a[5]; } 
  double& at31() { return a[6]; } 
  double& at32() { return a[7]; } 
  double& at33() { return a[8]; } 
  
  void rotate (Dvctr *v) const;
  void rotate (double *v) const;
  
  void copy_to_row_1 (const VectoR *v) { a[0] = v->x;  a[1] = v->y;  a[2] = v->z; }
  void copy_to_row_2 (const VectoR *v) { a[3] = v->x;  a[4] = v->y;  a[5] = v->z; }
  void copy_to_row_3 (const VectoR *v) { a[6] = v->x;  a[7] = v->y;  a[8] = v->z; }
  
  void copy_to_col_1 (const VectoR *v) { a[0] = v->x;  a[3] = v->y;  a[6] = v->z; }
  void copy_to_col_2 (const VectoR *v) { a[1] = v->x;  a[4] = v->y;  a[7] = v->z; }
  void copy_to_col_3 (const VectoR *v) { a[2] = v->x;  a[5] = v->y;  a[8] = v->z; }
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS ARRAY                          *** *** *** ***
//  ********************************************************************************************

#ifdef DEBUG
#define SRCCHCK    if (src==NULL) _errorr("source == NULL");
#define AVOID_CNST if (cnst)      _errorr("cnst")
#define AVOID_EXTA if (exta)      _errorr("exta")
#else
#define SRCCHCK
#define AVOID_CNST
#define AVOID_EXTA
#endif


enum arrayClassType { DvctrClass, LvctrClass, DscalClass };
enum arrayTypedef { ATlong , ATdouble };

class Dmtrx;

class Array
{
 public:
  Array () { }
  virtual ~Array () { }
  
  // give type def
  virtual arrayTypedef give_arrayTypedef (void) const = 0;
  
  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force=false) const = 0;
#endif
  
};



//  ********************************************************************************************
//  *** *** *** ***                        CLASS ARRAY1D                         *** *** *** ***
//  ********************************************************************************************
class Array1d : public Array
{
 public:
  Array1d () : Array() { }
  virtual ~Array1d() {  }
  
  virtual long give_size (void) const = 0;
  
  //  ***  PRINT  ***
  virtual int length_printed        (int precision) const = 0;
  virtual void print         (char *stream, int precision, double absZero=0.0) const = 0;
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS XSCAL                          *** *** *** ***
//  ********************************************************************************************
class Xscal : public Array1d
{
 public:
  Xscal () : Array1d() { }
  virtual ~Xscal() {  }
  
  long give_size (void) const { return 1; }
};

//  ********************************************************************************************
//  *** *** *** ***                         CLASS DSCAL                          *** *** *** ***
//  ********************************************************************************************
class Dscal : public Xscal
{
 public:
  double a;    // value
  
  Dscal (void)       : Xscal () { a = 0; }
  Dscal (double val) : Xscal () { a = val; }
  virtual ~Dscal() { }
  
  arrayClassType give_classid() const { return DscalClass; }
  arrayTypedef give_arrayTypedef (void) const { return ATdouble; }
  
  //  ***  OPERATORS  ***
  double  operator() (void) const { return a; }
  double& operator() (void)       { return a; }
  
  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force=false) const { fprintf (stdout, "Double scalar:  % .*e\n", PYP, a); }
#endif
  virtual int length_printed (int precision) const { return  1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT + 1; }
  virtual void print         (char *stream, int precision, double absZero=0.0) const;
};



//  ********************************************************************************************
//  *** *** *** ***                         CLASS XVCTR                          *** *** *** ***
//  ********************************************************************************************
class Xvctr : public Array1d
{
 protected:
  bool exta;    // external allocation, we cannot delete "a"
  bool cnst;    // constant, values can not be changed
  long shft;    // a is shifted

  long size;    // length of array
  long asize;   // length of the allocated array


 public:
  Xvctr (long s) : Array1d() { exta = cnst = false;  shft = 0;  size = asize = s; }
  virtual ~Xvctr() {  }

  virtual arrayClassType give_classid() const = 0;

  //  ***  SET  ***
  virtual void cpat (long i, const Xvctr *p, long j) = 0;
  virtual void shift (long val) = 0;
  void deshift (void) { this->shift(-shft); }
  virtual bool scan (const char *&src) = 0;

  //  ***  GET  ***
  virtual long give_size (void) const { return size; }

  //  ***  COMPUTE  ***
  virtual double give_lenght (void) const = 0;
  double give_zero (double abszero, double relzero) const;

  //  ***  PRINT  ***
  //virtual int length_printed        (int precision) const = 0;
  virtual int length_printed_vector (int precision) const = 0;
  virtual int length_printed_tensor (int precision) const = 0;
  //virtual void print           (char *stream, int precision, double absZero=0.0                     ) const = 0;
  virtual void print_vector    (char *stream, int precision, double absZero=0.0, bool zerorest=false) const = 0;
  virtual void print_symtensor (char *stream, int precision, double absZero=0.0, bool zerorest=false) const = 0;
  virtual void print_symtensor (FILE *stream, int precision, double absZero=0.0, bool zerorest=false) const = 0;
  //virtual void print_member          (char *stream, int precision, long i) const = 0;
};

//  ********************************************************************************************
//  *** *** *** ***                         CLASS LVCTR                          *** *** *** ***
//  ********************************************************************************************
class Lvctr : public Xvctr
{
  long *a;
  
 public:
  Lvctr (void)                    : Xvctr (0)       { a = NULL; }
  Lvctr (long s)                  : Xvctr (s)       { a = new long[asize]; }
  Lvctr (        const Lvctr *p)  : Xvctr (p->size) { a = new long[asize];  memcpy (a, p->a, size*sizeof(long)); }
  Lvctr (long s, const long   *p) : Xvctr (s)       { a = new long[asize];  memcpy (a, p,    size*sizeof(long)); }
  virtual ~Lvctr() { deshift(); if (a && !exta)  delete [] a; }
  
  arrayClassType give_classid() const { return LvctrClass; }
  arrayTypedef give_arrayTypedef (void) const { return ATlong; }
  
  //  ***  OPERATORS  ***
  long  operator[] (long i) const {
#ifdef DEBUG
    if (i<0 || i>=size)  _errorr ("i out of the range");
#endif
    return a[i];
  }
  
  long& operator[] (long i)       {
#ifdef DEBUG
    if (i<0 || i>=size)
      _errorr ("i out of the range");
#endif
    return a[i];
  }
  
  Lvctr& operator= (const Lvctr &src) {
    AVOID_CNST; AVOID_EXTA;
    if ( this == &src )  _errorr("oneself assignment");

    size = 0;
    this->realloc (src.size);
    size = src.size;
    
    memcpy (a, src.a, size*sizeof(long));
    
    return * this;
  }
  
  //  ***  SET  ***
  void add        (long val) { resize_preserve_vals(size+1);  a[size-1] = val; }
  void add_unique (long val) { if (! this->is_member(val))  this->add(val); }
  bool is_member (long val) const { return is_member_of_array (val, size, a); }
  Lvctr* free (void);
  void realloc (long newsize);
  Lvctr* resize_preserve_vals (long newsize) { AVOID_CNST; AVOID_EXTA;           this->realloc    (newsize);  size = newsize; return this; }
  Lvctr* resize_ignore_vals   (long newsize) { AVOID_CNST; AVOID_EXTA; size = 0; this->realloc    (newsize);  size = newsize; return this; }
  Lvctr* resize_to_asize      (void)         { AVOID_CNST; AVOID_EXTA;                                        size = asize;   return this; }
  Lvctr* assign_array (      long *array, long s);
  Lvctr* assign_array (const long *array, long s);
  void cpat (long i, const Xvctr *p, long j) {
    const Lvctr *pt = dynamic_cast<const Lvctr *>(p);
    if (!pt)  _errorr("");
    a[i] = pt->a[j];
  }
  void shift (long val) { AVOID_CNST; if (a) { shft += val;  a += val;  } }
  bool scan (const char *&src) { return SP_scan_array (src, size, a); }
  void fillYourselfBy (int val) { memset (a, val, size*sizeof(long)); }
  
  //  ***  GET  ***
  long*       give_ptr2val (long i=0)       { AVOID_CNST;  return a + i; }
  const long* give_ptr2val (long i=0) const {              return a + i; }
  
  
  //  ***  COMPUTE  ***
  double give_lenght (void) const;
  long give_sum                (void) const;
  long give_number_of_nonzeros (void) const;
  long give_number_of_zeros    (void) const;
  
  
  //  ***  EDIT  ***
  void zero (void) { AVOID_CNST;  memset (a,0,asize*sizeof(long)); }
  //
  //
  //
  //

  //
  //
  //
  //
  //

  //void addtms (const Dvctr &src, double tms);
  //void tmsby (double val);
  //void dvdby (double val);


  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force) const;
#endif
  //
  int  length_printed        (int precision) const;
  int  length_printed_vector (int precision) const;
  int  length_printed_tensor (int precision) const;
  void print           (char *stream, int precision, double absZero=0.0                     ) const;
  void print_vector    (char *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  void print_symtensor (char *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  void print_symtensor (FILE *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  //void print_member          (char *stream, int precision, long i) const { sprintf (stream, "%ld", a[i]); }
  //
  // void print_symtensor (FILE *stream) const;
  // void print_r2z_symtensor (FILE *stream, double relZero, double absZero) const;
  // void print_symtensor (char *stream) const;
  // void print_r2z_symtensor (char *stream, double relZero, double absZero) const;
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS DVCTR                          *** *** *** ***
//  ********************************************************************************************

class Dvctr : public Xvctr
{
 public:
  double *a;    // array of values

  Dvctr (void)                      : Xvctr (0)         {          a = NULL; }
  Dvctr (long s)                    : Xvctr (s)         {          a = new double[asize]; }
  Dvctr (        const Dvctr  *src) : Xvctr (src->size) { SRCCHCK; a = new double[asize];  memcpy (a, src->a, size*sizeof(double)); }
  Dvctr (long s, const double *src) : Xvctr (s)         { SRCCHCK; a = new double[asize];  memcpy (a, src,    size*sizeof(double)); }
  Dvctr (        const VectoR *src) : Xvctr (3)         { SRCCHCK; a = new double[asize];  a[0]=src->x;  a[1]=src->y;  a[2]=src->z; }
  virtual ~Dvctr() { deshift(); if (a && !exta)  delete [] a; }

  arrayClassType give_classid() const { return DvctrClass; }
  arrayTypedef give_arrayTypedef (void) const { return ATdouble; }
  
  //  ***  OPERATORS  ***
  double  operator[] (long i) const {
#ifdef DEBUG
    if (i<0 || i>=size)  _errorr ("i out of the range");
#endif
    return a[i];
  }
  double& operator[] (long i)       {
#ifdef DEBUG
    if (i<0 || i>=size)
      _errorr ("i out of the range");
#endif
    return a[i];
  }


  //  ***  SET  ***
  Dvctr* free (void);
  void realloc (long newsize);
  // resize, keep values
  Dvctr* resize_preserve_vals (long newsize) { AVOID_CNST; AVOID_EXTA;           this->realloc    (newsize);  size = newsize; return this; }
  // resize, ignore values
  Dvctr* resize_ignore_vals   (long newsize) { AVOID_CNST;
    //AVOID_EXTA;
    if (exta)
      _errorr("exta");
    size = 0; this->realloc    (newsize);  size = newsize; return this; }
  // resize, ignore values
  Dvctr* resize_to_asize      (void)         { AVOID_CNST; AVOID_EXTA;                                        size = asize;   return this; }
  //Dvctr* resizeup   (long newsize) { AVOID_CNST; AVOID_EXTA; this->reallocup  (newsize);  size = newsize; }
  //Dvctr* resizedown (long newsize) { AVOID_CNST; AVOID_EXTA; this->reallocdown(newsize);  size = newsize; }

  Dvctr* assign_array (      double *array, long s);
  Dvctr* assign_array (const double *array, long s);
  void cpat (long i, const Xvctr *p, long j) {
    const Dvctr *pt = dynamic_cast<const Dvctr *>(p);
    if (!pt)  _errorr("");
    a[i] = pt->a[j];
  }
  void shift (long val) { AVOID_CNST; if (a) { shft += val;  a += val;  } }
  bool scan (const char *&src) { return SP_scan_array (src, size, a); }

  //  ***  GET  ***
  double*       give_ptr2val (long i=0)       { AVOID_CNST;  return a + i; }
  const double* give_ptr2val (long i=0) const {              return a + i; }


  //  ***  COMPUTE  ***
  double give_lenght (void) const;
  double give_sum (void) const;
  long give_number_of_nonzeros (void) const;
  
  
  //  ***  EDIT  ***
  void zero (void) { AVOID_CNST;  memset (a,0,size*sizeof(double)); }
  void beCopyOf (const PoinT  &src)         { this->resize_ignore_vals(3);       a[0] = src.x;  a[1] = src.y;  a[2] = src.z; }
  void beCopyOf (const double *src)         { AVOID_CNST;                           memcpy (a, src,    size*sizeof(double)); }
  void beCopyOf (const double *src, long s) { this->resize_ignore_vals(s        );  memcpy (a, src,    size*sizeof(double)); }
  void beCopyOf (const Dvctr  &src)         { this->resize_ignore_vals(src. size);  memcpy (a, src.a,  size*sizeof(double)); }
  void beCopyOf (const Dvctr  *src)         { this->resize_ignore_vals(src->size);  memcpy (a, src->a, size*sizeof(double)); }
  
  void copy_to (double *dest) const { memcpy (dest, a, size*sizeof(double)); }
  
  void append (const Dvctr *src) {
    AVOID_CNST;
    this->resize_preserve_vals(size + src->size);
    memcpy (a+size-src->size, src->a, (size-src->size)*sizeof(double));
  }

  void be_tnsr (const Dmtrx &src);
  void be_vectproduct (const Dvctr &v1, const Dvctr &v2);
  void be_mean_of (const Dmtrx *src);
  double give_dotproduct (const Dvctr &v) const;
  
  
  void add    (const Dvctr &src);                    // this +=       src
  void addtms (const Dvctr &src, double tms);        // this += tms * src
  void sbt    (const Dvctr &src);                    // this -=       src
  void tmsby  (double val);
  void dvdby  (double val);

  void tnsrRotWith  (const Dmtrx &rot);
  void tnsrRotWith1 (const Dmtrx &rot);
  void tnsrRotWith2 (const Dmtrx &rot);
  void tnsrRotAxisZangle (double a);
  void normalize (void);
  double compute_squared_norm(void) const;
  double give_zero (double abszero, double relzero) const;

  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force) const;
#endif
  //
  int  length_printed        (int precision) const;
  int  length_printed_vector (int precision) const;
  int  length_printed_tensor (int precision) const;
  void print           (char *stream, int precision, double absZero=0.0                     ) const;
  void print_vector    (char *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  void print_symtensor (char *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  void print_symtensor (FILE *stream, int precision, double absZero=0.0, bool zerorest=false) const;
  
//  void print_vector          (char *stream, int precision) const;
//  void print_vector_or_zero  (char *stream, int precision) const;
//  void print_vector_zerorest (char *stream, int precision) const;
//  void print_member          (char *stream, int precision, long i) const { sprintf (stream, "% .*e", precision, a[i]); }
//  //
//  void print_symtensor      (FILE *stream, int precision) const;
//  void print_symtensor_zero (FILE *stream, int precision) const;
//  void print_r2z_symtensor  (FILE *stream, int precision, double relZero, double absZero) const;
//  void print_symtensor      (char *stream, int precision) const;
//  void print_symtensor_zero (char *stream, int precision) const;
//  void print_r2z_symtensor  (char *stream, int precision, double relZero, double absZero) const;
};





//  ********************************************************************************************
//  *** *** *** ***                         CLASS XMTRX                          *** *** *** ***
//  ********************************************************************************************
class Xmtrx : public Array
{
 protected:
  long crow;    // count of rows
  long ccol;    // count of collumns
  long asize;   // length of the allocated array

 public:
  Xmtrx (void)           : Array() {  crow = 0;  ccol = 0;  asize = crow*ccol; }
  Xmtrx (long r, long c) : Array() {  crow = r;  ccol = c;  asize = crow*ccol; }
  virtual ~Xmtrx() { }
  
  //  ***  SET  ***
  //  ***  GET  ***
  long give_crows (void) const { return crow; }
  long give_ccols (void) const { return ccol; }
};



//  ********************************************************************************************
//  *** *** *** ***                         CLASS LMTRX                          *** *** *** ***
//  ********************************************************************************************
class Lmtrx : public Xmtrx
{
 public:
  long *a;    // array of values

  Lmtrx (void)           : Xmtrx()    { a = NULL; }
  Lmtrx (long r, long c) : Xmtrx(r,c) { a = new long[asize]; }
  virtual ~Lmtrx() { if (a)  delete [] a; }
  
  arrayTypedef give_arrayTypedef (void) const { return ATlong; }
  
  //  ***  SET  ***
  void resize_ignore_vals (long r, long c);
  
  //  ***  GET  ***
  long       &operator()(long r, long c)       { return a[r*ccol + c]; }
  const long &operator()(long r, long c) const { return a[r*ccol + c]; }
  long*       give_ptr2val (long r=0, long c=0)       { return a + r*ccol + c; }
  const long* give_ptr2val (long r=0, long c=0) const { return a + r*ccol + c; }
  
  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force=false) const;
#endif
  
};


//  ********************************************************************************************
//  *** *** *** ***                         CLASS DMTRX                          *** *** *** ***
//  ********************************************************************************************
class Dmtrx : public Xmtrx
{
 protected:
  double *a;    // array of values
  
 public:
  Dmtrx (void)                            : Xmtrx()                    {          a = NULL; }
  Dmtrx (long r, long c)                  : Xmtrx(r,c)                 {          a = new double[asize]; }
  Dmtrx (               const Dmtrx *src) : Xmtrx(src->crow,src->ccol) { SRCCHCK; a = new double[asize];  memcpy (a, src->a, crow*ccol*sizeof(double)); }
  virtual ~Dmtrx() { delete [] a; }
  
  arrayTypedef give_arrayTypedef (void) const { return ATdouble; }
  
  //  ***  SET  ***
  void resize_ignore_vals (long r, long c);
  
  //  ***  GET  ***
  double       &operator()(long r, long c)       { return a[r*ccol + c]; }
  const double &operator()(long r, long c) const { return a[r*ccol + c]; }
  double*       give_ptr2val (long r=0, long c=0)       { return a + r*ccol + c; }
  const double* give_ptr2val (long r=0, long c=0) const { return a + r*ccol + c; }
  
  
  //  ***  EDIT  ***
  void zero (void) { memset (a, 0, crow*ccol*sizeof(double)); }
# ifndef DEBUG
  void copy_row (long r, const Dvctr &v) { if (ccol != v.give_size())  _errorr("size mismatch");  memcpy (a + r*ccol, v.give_ptr2val(), ccol*sizeof(double)); }
# else
  void copy_row (long r, const Dvctr &v) {                                                       memcpy (a + r*ccol, v.give_ptr2val(), ccol*sizeof(double)); }
# endif
  void copy_row (long r, const Elem3D *v) { a[r*ccol] = v->x;  a[r*ccol+1] = v->y;  a[r*ccol+2] = v->z; }
  void beCopyOf (const double *src)  { memcpy (a, src, crow*ccol*sizeof(double)); }
  void beCopyOf (const Dmtrx  &src)  { this->beCopyOf(src. give_ptr2val(0,0)); }
  void beCopyOf (const Dmtrx  *src)  { this->beCopyOf(src->give_ptr2val(0,0)); }
  void be_tnsr (const Dvctr &src);
  void be_transposition (const Dmtrx &src);
  void be_mtrxMmtrx (const Dmtrx &A, const Dmtrx &B);
  void tmsLeftBy (const Dmtrx &src);
  void tmsRightBy (const Dmtrx &src);
  
  //  ***  COMPUTE  ***
  double give_determinant (void) const;
  bool GaussSolve (Dvctr &x, const Dvctr &rhs);
  bool GaussSolve (Dvctr *x, const Dvctr *rhs, long nrhs);
  
  //  ***  PRINT  ***
#ifdef DEBUG
  virtual void printYourself (bool force=false) const;
#endif
  //void print_tensor     (FILE *stream, int precision) const;
  void print_tensor (FILE *stream, int precision, double absZero, bool zerorest) const;
  
};


} // namespace gelibspace

#endif // GELIB_ARRAYS_H