arrays.cpp

Go to the documentation of this file.

#include "arrays.h"

#include "gelib.h"

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


namespace gelibspace {

double give_copy_r2z (double a, double zero)
{
  return (-zero < a && a < zero) ? 0.0 : a ;
}
double* give_copy_r2z (double *a, long n, double zero)
{
  if (zero == 0.0)  return a;
  double *p = new double[n];
  while (n--)  p[n] = give_copy_r2z (a[n], zero);
  return p;
}


//  ********************************************************************************************
//  *** *** *** ***                         CLASS MatriX                         *** *** *** ***
//  ********************************************************************************************
void MatriX :: rotate (Dvctr *v) const
{
  long vs = v->give_size();
  if (vs != 3 && vs != 6) errol;
  
  ;            this->rotate(v->give_ptr2val(0));
  if (vs == 6) this->rotate(v->give_ptr2val(3));
}

void MatriX :: rotate (double *v) const
{
  double x = a[0] * v[0]  +  a[1] * v[1]  +  a[2] * v[2];
  double y = a[3] * v[0]  +  a[4] * v[1]  +  a[5] * v[2];
  double z = a[6] * v[0]  +  a[7] * v[1]  +  a[8] * v[2];
  
  v[0] = x;
  v[1] = y;
  v[2] = z;
}


//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS ELEM3D                         *** *** *** ***
//*  ********************************************************************************************

//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS POINT                          *** *** *** ***
//*  ********************************************************************************************
void PoinT :: beRotatedPoint (const Dmtrx *rot, const PoinT *point)
{
#ifdef DEBUG
  if (rot->give_crows() != 3 || rot->give_ccols() != 3)  errol;
#endif
  const double *r = rot->give_ptr2val();

  x = point->x * r[0] + point->y * r[1] + point->z * r[2];
  y = point->x * r[3] + point->y * r[4] + point->z * r[5];
  z = point->x * r[6] + point->y * r[7] + point->z * r[8];
}
void PoinT :: beRotatedBy (const double *r)
{
  double xx = x * r[0] + y * r[1] + z * r[2];
  double yy = x * r[3] + y * r[4] + z * r[5];
  double zz = x * r[6] + y * r[7] + z * r[8];
  x = xx;
  y = yy;
  z = zz;
}
void PoinT :: beRotatedBy2d (const double *r)
{
  double xx = x * r[0] + y * r[1];
  double yy = x * r[2] + y * r[3];
  x = xx;
  y = yy;
  z = 0.0;
}

//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS DSCAL                          *** *** *** ***
//*  ********************************************************************************************
void Dscal :: print (char *stream, int precision, double absZero) const
{
  sprintf (stream, "% .*e", precision, (absZero != 0.0) ? give_copy_r2z (a, absZero) : a);
  stream += 2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT;
}



double Xvctr :: give_zero (double abszero, double relzero) const
{
  relzero *= this->give_lenght();
  return (relzero > abszero ? relzero : abszero);
}

//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS LVCTR                          *** *** *** ***
//*  ********************************************************************************************

//*  ***  SET  ***
void Lvctr :: realloc (long newsize)
{
  if (newsize <= asize)  return;
  
  asize = newsize;
  long* newa = new long[asize];
  
  for (newsize=0; newsize<size; newsize++)  *newa++ = *a++;
  
  delete [] (a-size);
  a = newa-size;
}

Lvctr* Lvctr :: assign_array (long *array, long s)
{
  if (a)  errol;
  cnst = false;
  exta = true;
  size = s;
  asize = 0;
  a = array;
  return this;
}
Lvctr* Lvctr :: assign_array (const long *array, long s)
{
  if (a)  errol;
  cnst = true;
  exta = true;
  size = s;
  asize = 0;
  a = (long *)array;
  return this;
}
Lvctr* Lvctr :: free (void) {
  if (a && !exta)  delete [] a;
  a = NULL;
  exta = cnst = false;
  shft = size = asize = 0;
  return this;
}

//*  ***  COMPUTE  ***
double Lvctr :: give_lenght (void) const
{
  double length = 0.0;
  for (long i=0; i<size; i++)
    length += a[i] * a[i];
  
  return sqrt(length);
}
long Lvctr :: give_sum                (void) const { long sum = 0;  for (long i=0; i<size; i++)              sum += a[i];  return sum; }
long Lvctr :: give_number_of_nonzeros (void) const { long sum = 0;  for (long i=0; i<size; i++)  if ( a[i])  sum++;        return sum; }
long Lvctr :: give_number_of_zeros    (void) const { long sum = 0;  for (long i=0; i<size; i++)  if (!a[i])  sum++;        return sum; }


//*  ***
//*  ***  PRINT  ***
//*  ***
int Lvctr :: length_printed        (int precision) const { return size * (25+1); }
int Lvctr :: length_printed_vector (int precision) const { return    3 * (25+1); }
int Lvctr :: length_printed_tensor (int precision) const { return    9 * (25+1); }

//*  ***  VECTOR
#ifdef DEBUG
void Lvctr :: printYourself (bool force) const
{
  if (size > 20 && !force)  _errorr("Count of rows is > 20, to force printing set first parameter as true");
  
  fprintf (stdout, "Long vector [%ld]:\n", size);
  
  for (long i=0; i<size; i++)  fprintf (stdout, "%ld ", a[i]);
  
  fprintf (stdout, "\n");
}
#endif

void Lvctr :: print (char *stream, int precision, double absZero) const
{
  if (size < 1)  errol;
  
  int numlen = precision+1;
  
  ;                             sprintf (stream,  "%ld", a[0]);  stream += numlen-1;
  for (long i=1; i<size; i++) { sprintf (stream, " %*ld", precision, a[i]);  stream += numlen;  }
}
void Lvctr :: print_vector (char *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size < 0)  errol; }
  else if (size < 3)  errol;
  
  sprintf (stream, "%*ld %*ld %*ld",
       precision, size>=1 ? a[0] : 0,
       precision, size>=2 ? a[1] : 0,
       precision, size>=3 ? a[2] : 0);
}

//*  ***  TENSOR
void pprint_symtensor (char *stream, int precision, const long *a)
{
  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[0], precision, a[5], precision, a[4]);   while (stream[0] != '\0') stream++;
  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[5], precision, a[1], precision, a[3]);   while (stream[0] != '\0') stream++;
  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[4], precision, a[3], precision, a[2]);
}
void pprint_symtensor_zeroL (char *stream, int precision)
{
  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);   while (stream[0] != '\0') stream++;
  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);   while (stream[0] != '\0') stream++;
  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
}
void pprint_symtensor (FILE *stream, int precision, const long *a)
{
  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[0], precision, a[5], precision, a[4]);
  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[5], precision, a[1], precision, a[3]);
  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[4], precision, a[3], precision, a[2]);
}
void pprint_symtensor_zeroL (FILE *stream, int precision)
{
  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
}

void Lvctr :: print_symtensor (char *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size != 0  &&  size != 6)  errol; }
  else if (size != 6)  errol;
  
  if (size) pprint_symtensor      (stream, precision, a);
  else      pprint_symtensor_zeroL (stream, precision);
}
void Lvctr :: print_symtensor (FILE *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size != 0  &&  size != 6)  errol; }
  else if (size != 6)  errol;
  
  if (size) pprint_symtensor      (stream, precision, a);
  else      pprint_symtensor_zeroL (stream, precision);
}



//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS DVCTR                          *** *** *** ***
//*  ********************************************************************************************
//*  ***  SET  ***
void Dvctr :: realloc (long newsize)
{
  if (newsize <= asize)  return;
  
  asize = newsize;
  double* newa = new double[asize];
  
  for (newsize=0; newsize<size; newsize++)  *newa++ = *a++;
  
  delete [] (a-size);
  a = newa-size;
}

Dvctr* Dvctr :: assign_array (double *array, long s)
{
  if (a)  errol;
  cnst = false;
  exta = true;
  size = s;
  asize = 0;
  a = array;
  return this;
}
Dvctr* Dvctr :: assign_array (const double *array, long s)
{
  if (a)
    errol;
  cnst = true;
  exta = true;
  size = s;
  asize = 0;
  a = (double *)array;
  return this;
}
Dvctr* Dvctr :: free (void) {
  if (a && !exta)  delete [] a;
  a = NULL;
  exta = cnst = false;
  shft = size = asize = 0;
  return this;
}

//*  ***  COMPUTE  ***
double Dvctr :: give_lenght (void) const
{
# ifdef DEBUG
  if (size == 0)  _warningg ("size = 0 in give_length function");
# endif
  
  double length = 0.0;
  for (long i=0; i<size; i++)
    length += a[i] * a[i];
  
  return sqrt(length);
}
double Dvctr :: give_sum (void) const
{
  double sum = 0;
  for (long i=0; i<size; i++)
    sum += a[i];
  
  return sum;
}

//*  ***  EDIT  ***
void Dvctr :: be_tnsr (const Dmtrx &src)
{
  AVOID_CNST;
  if ( src.give_crows() != 3  &&  src.give_ccols() != 3 )  _errorr ("size of mtrx != 3x3");
  
  const double *asrc = src.give_ptr2val(0,0);
  this->resize_ignore_vals(6);

  a[0] = asrc[0];
  a[1] = asrc[4];
  a[2] = asrc[8];
  
  a[3] = asrc[5];
  a[4] = asrc[2];
  a[5] = asrc[1];
}

void Dvctr :: be_vectproduct (const Dvctr &v1, const Dvctr &v2)
{
  AVOID_CNST;
  if ( v1.give_size() != 3  ||  v2.give_size() != 3 )
    _errorr ("size of vectors is not equal to 3");
  
  this->resize_ignore_vals(3);
  
  const double *a1 = v1.give_ptr2val();
  const double *a2 = v2.give_ptr2val();
  
  a[0] = a1[1] * a2[2] - a2[1] * a1[2];
  a[1] = a1[2] * a2[0] - a2[2] * a1[0];
  a[2] = a1[0] * a2[1] - a2[0] * a1[1];
}

void Dvctr :: be_mean_of (const Dmtrx *src)
{
  long    n    = src->give_crows();
  long    m    = src->give_ccols();
  const double *psrc = src->give_ptr2val();
  
  this->beCopyOf(psrc, m);
  
  if (n>1) {
    psrc += m;
    while (--n)
      for (long i=0; i<m; i++)
    a[i] += *(psrc++);
    
    this->dvdby(src->give_crows());
  }
}

double Dvctr :: give_dotproduct (const Dvctr &v) const
{
  double prd = 0.0;
  for (long i=0; i<this->size; i++)
    prd += a[i] * v.a[i];
  
  return prd;
}
  

void Dvctr :: add    (const Dvctr &src)             { AVOID_CNST;  if (size != src.give_size()) _errorr("size mismatch");  const double *asrc = src.give_ptr2val();  for (long i=0; i<size; i++)  a[i] +=       asrc[i]; }
void Dvctr :: addtms (const Dvctr &src, double tms) { AVOID_CNST;  if (size != src.give_size()) _errorr("size mismatch");  const double *asrc = src.give_ptr2val();  for (long i=0; i<size; i++)  a[i] += tms * asrc[i]; }
void Dvctr :: sbt    (const Dvctr &src)             { AVOID_CNST;  if (size != src.give_size()) _errorr("size mismatch");  const double *asrc = src.give_ptr2val();  for (long i=0; i<size; i++)  a[i] -=       asrc[i]; }
void Dvctr :: dvdby (double val)
{
  AVOID_CNST;
  for (long i=0; i<size; i++)
    a[i] = a[i] / val;
}
void Dvctr :: tmsby (double val)
{
  AVOID_CNST;
  for (long i=0; i<size; i++)
    a[i] = a[i] * val;
}
//void Dvctr :: tnsrRotWith1 (const Dmtrx &rot)
//{
//  if ( size != 6 )  _errorr ("size of vectors is not equal to 6");
//  
//  Dmtrx sss, rotT;
//  
//  sss.be_tnsr(*this);
//  rotT.be_transposition(rot);
//  
//  sss.tmsLeftBy  (rot);
//  sss.tmsRightBy (rotT);
//  
//  this->be_tnsr(sss);
//}
void Dvctr :: tnsrRotWith (const Dmtrx &rot)
{
  AVOID_CNST;
  if ( size != 6 )  _errorr ("size of vectors is not equal to 6");
  
  double aa[6];
  
  aa[0] = rot(0,0)*rot(0,0) * a[0]  +  rot(0,1)*rot(0,1) * a[1]  +  rot(0,2)*rot(0,2) * a[2]  +  2*rot(0,0)*rot(0,1) * a[5]  +  2*rot(0,0)*rot(0,2) * a[4]  +  2*rot(0,1)*rot(0,2) * a[3];
  aa[1] = rot(1,0)*rot(1,0) * a[0]  +  rot(1,1)*rot(1,1) * a[1]  +  rot(1,2)*rot(1,2) * a[2]  +  2*rot(1,0)*rot(1,1) * a[5]  +  2*rot(1,0)*rot(1,2) * a[4]  +  2*rot(1,1)*rot(1,2) * a[3];
  aa[2] = rot(2,0)*rot(2,0) * a[0]  +  rot(2,1)*rot(2,1) * a[1]  +  rot(2,2)*rot(2,2) * a[2]  +  2*rot(2,0)*rot(2,1) * a[5]  +  2*rot(2,0)*rot(2,2) * a[4]  +  2*rot(2,1)*rot(2,2) * a[3];
  
  aa[5] = rot(0,0)*rot(1,0) * a[0]  +  rot(0,1)*rot(1,1) * a[1]  +  rot(0,2)*rot(1,2) * a[2]  +  (rot(0,0)*rot(1,1) + rot(0,1)*rot(1,0)) * a[5]  +  (rot(0,1)*rot(1,2) + rot(0,2)*rot(1,1)) * a[3]  +  (rot(0,0)*rot(1,2) + rot(0,2)*rot(1,0)) * a[4];
  aa[3] = rot(1,0)*rot(2,0) * a[0]  +  rot(1,1)*rot(2,1) * a[1]  +  rot(1,2)*rot(2,2) * a[2]  +  (rot(1,0)*rot(2,1) + rot(1,1)*rot(2,0)) * a[5]  +  (rot(1,1)*rot(2,2) + rot(1,2)*rot(2,1)) * a[3]  +  (rot(1,0)*rot(2,2) + rot(1,2)*rot(2,0)) * a[4];
  aa[4] = rot(0,0)*rot(2,0) * a[0]  +  rot(0,1)*rot(2,1) * a[1]  +  rot(0,2)*rot(2,2) * a[2]  +  (rot(0,0)*rot(2,1) + rot(0,1)*rot(2,0)) * a[5]  +  (rot(0,1)*rot(2,2) + rot(0,2)*rot(2,1)) * a[3]  +  (rot(0,0)*rot(2,2) + rot(0,2)*rot(2,0)) * a[4];
  
  this->beCopyOf(aa);
}

//void Dvctr :: tnsrRotAxisZangle (double a)
//{
//  Dmtrx rot(3,3);
//  rot(0,0) = cos(a);  rot(0,1) = -sin(a);  rot(0,2) = 0;
//  rot(1,0) = sin(a);  rot(1,1) =  cos(a);  rot(1,2) = 0;
//  rot(2,0) =      0;  rot(2,1) =       0;  rot(2,2) = 1;
//  
//  this->tnsrRotWith(rot);
//}
void Dvctr :: tnsrRotAxisZangle (double alpha)
{
  AVOID_CNST;
  if ( size != 6 )  _errorr ("size of vectors is not equal to 6");
  
  double c = cos(alpha);
  double s = sin(alpha);
  
  double aa[6];
  
  aa[0] = c*c * a[0]  +  s*s * a[1]  -  2*c*s * a[5] ;
  aa[1] = s*s * a[0]  +  c*c * a[1]  +  2*s*c * a[5] ;
  aa[2] = a[2] ;
  
  aa[5] =  c*s * a[0]  - s*c * a[1]  +  (c*c - s*s) * a[5] ;
  aa[3] =  c   * a[3]  + s   * a[4];
  aa[4] =   -s * a[3]  +   c * a[4];
  
  this->beCopyOf(aa);
}



void Dvctr :: normalize (void)
{
  AVOID_CNST;
  double length = this->give_lenght();
  this->dvdby(length);
}

double Dvctr :: compute_squared_norm (void) const
{
  double *p = this->a;
  int     i = this->size;
  double norm2 = 0.0;
  while (i--) {
    norm2 += (*p) * (*p);
    p++;
  }
  
  return norm2;
}

//void Dvctr :: round2zero (double relZero, double absZero)
//{
//  AVOID_CNST;
//  double l;
//  l = this->give_lenght();
//  for (long i=0; i<size; i++)
//    if (fabs(a[i]/l) < relZero  ||  fabs(a[i]) < absZero)  a[i] = 0.0;
//  
//}


//void Dvctr :: round2zero (double relZero, double absZero)
//{
//  AVOID_CNST;
//  double l;
//  l = this->give_lenght();
//  for (long i=0; i<size; i++)
//    if (fabs(a[i]/l) < relZero  ||  fabs(a[i]) < absZero)  a[i] = 0.0;
//  
//}

//if (zero < 0.0) {
//    double l = 0.0;
//    for (long i=0; i<n; i++)  l += a[i] * a[i];
//    l = sqrt(l) * zero;
//    
//    if ()
//  




//*  ***
//*  ***   PRINT   ***
//*  ***
int Dvctr :: length_printed        (int precision) const { return size * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }
int Dvctr :: length_printed_vector (int precision) const { return    3 * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }
int Dvctr :: length_printed_tensor (int precision) const { return    9 * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }

//*  ***  VECTOR
#ifdef DEBUG
void Dvctr :: printYourself (bool force) const
{
  if (size > 20 && !force)  _errorr("Count of rows is > 20, to force printing set first parameter as true");
  
  fprintf (stdout, "Double vector [%ld]:\n", size);
  
  for (long i=0; i<size; i++)  fprintf (stdout, "% .*e ", PYP, a[i]);
  
  fprintf (stdout, "\n");
}
#endif

void Dvctr :: print (char *stream, int precision, double absZero) const
{
  if (size < 1)  errol;
  
  double *p = a;
  if (absZero != 0.0)  p = give_copy_r2z (a, size, absZero);
  
  int numlen = 1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT;
  
  ;                             sprintf (stream,  "% .*e", precision, p[0]);  stream += numlen-1;
  for (long i=1; i<size; i++) { sprintf (stream, " % .*e", precision, p[i]);  stream += numlen;  }
  
  if (p != a)  delete [] p;
}
void Dvctr :: print_vector (char *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size < 0)  errol; }
  else if (size < 3)  errol;
  
  double *p = a;
  if (absZero != 0.0 && size)  p = give_copy_r2z (a, 3, absZero);
  
  sprintf (stream, "% .*e % .*e % .*e",
       precision, size>=1 ? p[0] : 0.0,
       precision, size>=2 ? p[1] : 0.0,
       precision, size>=3 ? p[2] : 0.0);
  
  if (p != a)  delete [] p;
}

//*  ***  TENSOR

void pprint_tensor (FILE *stream, int precision, const double *a)
{
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[1], precision, a[2]);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[3], precision, a[4], precision, a[5]);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[6], precision, a[7], precision, a[8]);
}
void pprint_symtensor (char *stream, int precision, const double *a)
{
  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);   while (stream[0] != '\0') stream++;
  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);   while (stream[0] != '\0') stream++;
  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
}
void pprint_symtensor (FILE *stream, int precision, const double *a)
{
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
}

void pprint_tensor_zero (char *stream, int precision)
{
  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);   while (stream[0] != '\0') stream++;
  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);   while (stream[0] != '\0') stream++;
  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
}
void pprint_tensor_zero (FILE *stream, int precision)
{
  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
}



void Dvctr :: print_symtensor (char *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size != 0  &&  size != 6)  errol; }
  else if (size != 6)  errol;
  
  double *p = a;
  if (absZero != 0.0 && size)  p = give_copy_r2z (a, 6, absZero);
  
  if (size) pprint_symtensor   (stream, precision, p);
  else      pprint_tensor_zero (stream, precision);
  
  if (p != a)  delete [] p;
}
void Dvctr :: print_symtensor (FILE *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (size != 0  &&  size != 6)  errol; }
  else if (size != 6)  errol;
  
  double *p = a;
  if (absZero != 0.0 && size)  p = give_copy_r2z (a, 6, absZero);
  
  if (size) pprint_symtensor   (stream, precision, p);
  else      pprint_tensor_zero (stream, precision);
  
  if (p != a)  delete [] p;
}


//void Dvctr :: print_symtensor (char *stream, int precision) const
//{
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);   while (stream[0] != '\0') stream++;
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);   while (stream[0] != '\0') stream++;
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
//}
//void Dvctr :: print_symtensor_zero (char *stream, int precision) const
//{
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);   while (stream[0] != '\0') stream++;
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);   while (stream[0] != '\0') stream++;
//  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
//}
//void Dvctr :: print_symtensor (FILE *stream, int precision) const
//{
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
//}
//void Dvctr :: print_symtensor_zero (FILE *stream, int precision) const
//{
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
//}





//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS LMTRX                          *** *** *** ***
//*  ********************************************************************************************

//*  *************
//*  ***  SET  ***
//*  *************
void Lmtrx :: resize_ignore_vals (long r, long c)
{
  if (crow == r  &&  ccol == c)  return;
  
  crow = r;  ccol = c;
  
  if (asize < crow*ccol) {
    if (a)  delete [] a;
    asize = crow*ccol;
    a = new long[asize];
  }
}


//*  ***************
//*  ***  PRINT  ***
//*  ***************
#ifdef DEBUG
void Lmtrx :: printYourself (bool force) const
{
  if (crow > 50)  _warningg("count of rows     is > 50");
  if (ccol > 50)  _warningg("count of collumns is > 50");
  
  fprintf (stdout, "Long matrix [%ld,%ld]:\n", crow, ccol);
  
  long i, j;
  for (i=0; i<crow; i++) {
    for (j=0; j<ccol; j++)
      fprintf (stdout, "%ld", a[i*ccol + j]);
    
    fprintf (stdout, "\n");
  }
}
#endif



//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS DMTRX                          *** *** *** ***
//*  ********************************************************************************************

//*  *************
//*  ***  SET  ***
//*  *************
void Dmtrx :: resize_ignore_vals (long r, long c)
{
  if (crow == r  &&  ccol == c)  return;
  
  crow = r;  ccol = c;
  
  if (asize < crow*ccol) {
    if (a)  delete [] a;
    asize = crow*ccol;
    a = new double[asize];
  }
}


//*  **************
//*  ***  EDIT  ***
//*  **************
void Dmtrx :: be_tnsr (const Dvctr &src)
{
  if ( src.give_size() != 6 )  _errorr ("size of vectors is not equal to 6");
  
  const double *asrc = src.give_ptr2val();
  this->resize_ignore_vals(3,3);
  
  a[0]        = asrc[0];
  a[4]        = asrc[1];
  a[8]        = asrc[2];
  
  a[5] = a[7] = asrc[3];
  a[2] = a[6] = asrc[4];
  a[1] = a[3] = asrc[5];
}
void Dmtrx :: be_transposition (const Dmtrx &src)
{
  this->resize_ignore_vals(src.give_ccols(), src.give_crows());
  
  const double *asrc = src.give_ptr2val(0,0);
  long i, j;
  for (i=0; i<crow; i++)
    for (j=0; j<ccol; j++)
      a[i*ccol + j] = asrc[j*ccol + i];
  
}
void Dmtrx :: be_mtrxMmtrx (const Dmtrx &A, const Dmtrx &B)
{
  if ( A.give_ccols() != B.give_crows() )  _errorr("size mismatch");
  long ccr = A.give_ccols();
  
  this->resize_ignore_vals(A.give_crows(), B.give_ccols());
  this->zero();
  
  const double *aA = A.give_ptr2val(0,0);
  const double *aB = B.give_ptr2val(0,0);
  
  long i, j, k;
  for (i=0; i<crow; i++)
    for (j=0; j<ccol; j++)
      for (k=0; k<ccr; k++)
    a[i*ccol + j] += aA[i*ccr+k] * aB[k*ccol+j];
  
}
void Dmtrx :: tmsLeftBy (const Dmtrx &src)
{
  Dmtrx mtrxXthis;
  mtrxXthis.be_mtrxMmtrx (src, *this);
  
  this->beCopyOf(mtrxXthis);
}
void Dmtrx :: tmsRightBy (const Dmtrx &src)
{
  Dmtrx thisXmtrx;
  thisXmtrx.be_mtrxMmtrx (*this, src);
  
  this->beCopyOf(thisXmtrx);
}


//*  *****************
//*  ***  COMPUTE  ***
//*  *****************
double Dmtrx :: give_determinant (void) const
{
  if      (crow == ccol && ccol == 2)  return a[0] * a[3] - a[1] * a[2];
  else if (crow == ccol && ccol == 3)  return a[0] * a[4] * a[8] + a[3] * a[7] * a[2] +
                          a[6] * a[1] * a[5] - a[6] * a[4] * a[2] -
                              a[7] * a[5] * a[0] - a[8] * a[3] * a[1] ;
  else  _errorr3("cannot compute determinant of matrix %ld x %ld", crow, ccol);
  return 0.0;
}

bool Dmtrx :: GaussSolve (Dvctr &x, const Dvctr &rhs)
{
#ifdef DEBUG
  if (crow != ccol)             _errorr3 ("matrix is not square, %ld != %ld", crow, ccol);
  if (crow != rhs.give_size())  _errorr3 ("matrix and rhs dimension differs, %ld != %ld", crow, rhs.give_size());
#endif
  
  int i, j, k, maxrow;
  double tmp;
  
  x = rhs;
  
  for (i=0; i<crow-1; i++) {
    
    // Find the row with the largest first value
    maxrow = i;
    for (j=i+1; j<crow; j++)
      if (fabs(a[j*ccol+i]) > fabs(a[maxrow*ccol+i]))
    maxrow = j;
    
    // Singular matrix?
    if (fabs(a[maxrow*ccol+i]) < 1.0e-15)
      return(false);
    
    // Swap the maxrow and ith row
    if (maxrow != i) {
      for (j=i; j<crow; j++) {
    tmp = a[i*ccol+j];
    a[i*ccol+j] = a[maxrow*ccol+j];
    a[maxrow*ccol+j] = tmp;
      }
      
      tmp = x[i];
      x[i] = x[maxrow];
      x[maxrow] = tmp;
    }
    
    // Eliminate the ith element of the jth row
    for (j=i+1; j<crow; j++) {
      tmp = a[j*ccol+i] / a[i*ccol+i];
      
      for (k=i; k<crow; k++)
    a[j*ccol+k] -= a[i*ccol+k] * tmp;
      
      x[j] -= x[i] * tmp;
    }
  }
  
  // Singular matrix?
  if (fabs(a[crow*ccol-1]) < 1.0e-15)
    return(false);
  
  // Do the back substitution
  for (i=crow-1; i>=0; i--) {
    tmp = 0;
    for (j=i+1; j<crow; j++)
      tmp += a[i*ccol+j] * x[j];
    
    x[i] = (x[i] - tmp) / a[i*ccol+i];
  }
  
  return true;
}


bool Dmtrx :: GaussSolve (Dvctr *x, const Dvctr *rhs, long nrhs)
{
#ifdef DEBUG
  if (crow != ccol)                _errorr3 ("matrix is not square, %ld != %ld", crow, ccol);
  if (crow != rhs[0].give_size())  _errorr3 ("matrix and rhs dimension differs, %ld != %ld", crow, rhs[0].give_size());
#endif
  
  int i, j, k, n, maxrow;
  double tmp;
  
  for (n=0; n<nrhs; n++)  x[n].beCopyOf(rhs[n]);
  
  for (i=0; i<crow-1; i++) {
    
    // Find the row with the largest first value
    maxrow = i;
    for (j=i+1; j<crow; j++)
      if (fabs(a[j*ccol+i]) > fabs(a[maxrow*ccol+i]))
    maxrow = j;
    
    // Singular matrix?
    if (fabs(a[maxrow*ccol+i]) < 1.0e-15)
      return(false);
    
    // Swap the maxrow and ith row
    if (maxrow != i) {
      for (j=i; j<crow; j++) {
    tmp = a[i*ccol+j];
    a[i*ccol+j] = a[maxrow*ccol+j];
    a[maxrow*ccol+j] = tmp;
      }
      
      for (n=0; n<nrhs; n++) {
    tmp = x[n][i];
    x[n][i] = x[n][maxrow];
    x[n][maxrow] = tmp;
      }
    }
    
    // Eliminate the ith element of the jth row
    for (j=i+1; j<crow; j++) {
      tmp = a[j*ccol+i] / a[i*ccol+i];
      
      for (k=i; k<crow; k++)
    a[j*ccol+k] -= a[i*ccol+k] * tmp;
      
      for (n=0; n<nrhs; n++)
    x[n][j] -= x[n][i] * tmp;
    }
  }
  
  // Singular matrix?
  if (fabs(a[crow*ccol-1]) < 1.0e-15)
    return(false);
  
  // Do the back substitution
  for (i=crow-1; i>=0; i--) {
    for (n=0; n<nrhs; n++) {
      tmp = 0;
      for (j=i+1; j<crow; j++)
    tmp += a[i*ccol+j] * x[n][j];
      
      x[n][i] = (x[n][i] - tmp) / a[i*ccol+i];
    }
  }
  
  return true;
}


//*  ***************
//*  ***  PRINT  ***
//*  ***************
#ifdef DEBUG
void Dmtrx :: printYourself (bool force) const
{
  if (crow > 10)  _errorr("count of rows is > 10");
  if (ccol > 10)  _errorr("count of collumns is > 10");
  
  fprintf (stdout, "Double matrix [%ld,%ld]:\n", crow, ccol);
  
  long i, j;
  for (i=0; i<crow; i++) {
    for (j=0; j<ccol; j++)
      fprintf (stdout, "% .*e ", PYP, a[i*ccol + j]);
    
    fprintf (stdout, "\n");
  }
}
#endif


//void Dmtrx :: print_tensor (FILE *stream, int precision) const
//{
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[1], precision, a[2]);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[3], precision, a[4], precision, a[5]);
//  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[6], precision, a[7], precision, a[8]);
//}

void Dmtrx :: print_tensor (FILE *stream, int precision, double absZero, bool zerorest) const
{
  if (zerorest) { if (! ((crow == 0 && ccol == 0)  ||  (crow == 3 && ccol == 3)) )  errol; }
  else if (crow != 3 || ccol != 3)  errol;
  
  double *p = a;
  if (absZero != 0.0 && crow)  p = give_copy_r2z (a, 9, absZero);
  
  if (crow) pprint_tensor      (stream, precision, p);
  else      pprint_tensor_zero (stream, precision);
  
  if (p != a)  delete [] p;
}


} // namespace gelibspace