MIDAS  0.75
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
arrays.cpp
Go to the documentation of this file.
1 #include "arrays.h"
2 
3 #include "gelib.h"
4 
5 #include <stdlib.h>
6 #include <string.h>
7 #include <math.h>
8 
9 
10 namespace midaspace {
11 
13 double give_copy_r2z (double a, double zero)
14 {
15  return (-zero < a && a < zero) ? 0.0 : a ;
16 }
17 double* give_copy_r2z (double *a, long n, double zero)
18 {
19  if (zero == 0.0) return a;
20  double *p = new double[n];
21  while (n--) p[n] = give_copy_r2z (a[n], zero);
22  return p;
23 }
24 
25 
26 // ********************************************************************************************
27 // *** *** *** *** CLASS MatriX *** *** *** ***
28 // ********************************************************************************************
30 void MatriX :: rotate (Dvctr *v) const
31 {
32  long vs = v->give_size();
33  if (vs != 3 && vs != 6) errol;
34 
35  ; this->rotate(v->give_ptr2val(0));
36  if (vs == 6) this->rotate(v->give_ptr2val(3));
37 }
38 
40 void MatriX :: rotate (double *v) const
41 {
42  double x = a[0] * v[0] + a[1] * v[1] + a[2] * v[2];
43  double y = a[3] * v[0] + a[4] * v[1] + a[5] * v[2];
44  double z = a[6] * v[0] + a[7] * v[1] + a[8] * v[2];
45 
46  v[0] = x;
47  v[1] = y;
48  v[2] = z;
49 }
50 
51 
52 //* ********************************************************************************************
53 //* *** *** *** *** CLASS ELEM3D *** *** *** ***
54 //* ********************************************************************************************
55 
56 //* ********************************************************************************************
57 //* *** *** *** *** CLASS POINT *** *** *** ***
58 //* ********************************************************************************************
60 void PoinT :: beRotatedPoint (const Dmtrx *rot, const PoinT *point)
61 {
62 #ifdef DEBUG
63  if (rot->give_crows() != 3 || rot->give_ccols() != 3) errol;
64 #endif
65  const double *r = rot->give_ptr2val();
66 
67  x = point->x * r[0] + point->y * r[1] + point->z * r[2];
68  y = point->x * r[3] + point->y * r[4] + point->z * r[5];
69  z = point->x * r[6] + point->y * r[7] + point->z * r[8];
70 }
71 
72 
73 
74 //* ********************************************************************************************
75 //* *** *** *** *** CLASS DSCAL *** *** *** ***
76 //* ********************************************************************************************
77 void Dscal :: print (char *stream, int precision, double absZero) const
78 {
79  sprintf (stream, "% .*e", precision, (absZero != 0.0) ? give_copy_r2z (a, absZero) : a);
80  stream += 2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT;
81 }
82 
83 
84 
85 double Xvctr :: give_zero (double abszero, double relzero) const
86 {
87  relzero *= this->give_lenght();
88  return (relzero > abszero ? relzero : abszero);
89 }
90 
91 //* ********************************************************************************************
92 //* *** *** *** *** CLASS LVCTR *** *** *** ***
93 //* ********************************************************************************************
94 
95 //* *** SET ***
97 void Lvctr :: realloc (long newsize)
98 {
99  if (newsize <= asize) return;
100 
101  asize = newsize;
102  long* newa = new long[asize];
103 
104  for (newsize=0; newsize<size; newsize++) *newa++ = *a++;
105 
106  delete [] (a-size);
107  a = newa-size;
108 }
109 
111 Lvctr* Lvctr :: assign_array (long *array, long s)
112 {
113  if (a) errol;
114  cnst = false;
115  exta = true;
116  size = s;
117  asize = 0;
118  a = array;
119  return this;
120 }
121 Lvctr* Lvctr :: assign_array (const long *array, long s)
122 {
123  if (a) errol;
124  cnst = true;
125  exta = true;
126  size = s;
127  asize = 0;
128  a = (long *)array;
129  return this;
130 }
131 Lvctr* Lvctr :: free (void) {
132  if (a && !exta) delete [] a;
133  a = NULL;
134  exta = cnst = false;
135  shft = size = asize = 0;
136  return this;
137 }
138 
139 //* *** COMPUTE ***
140 double Lvctr :: give_lenght (void) const
141 {
142  double length = 0.0;
143  for (long i=0; i<size; i++)
144  length += a[i] * a[i];
145 
146  return sqrt(length);
147 }
148 long Lvctr :: give_sum (void) const { long sum = 0; for (long i=0; i<size; i++) sum += a[i]; return sum; }
149 long Lvctr :: give_number_of_nonzeros (void) const { long sum = 0; for (long i=0; i<size; i++) if ( a[i]) sum++; return sum; }
150 long Lvctr :: give_number_of_zeros (void) const { long sum = 0; for (long i=0; i<size; i++) if (!a[i]) sum++; return sum; }
151 
152 
153 //* ***
154 //* *** PRINT ***
155 //* ***
156 int Lvctr :: length_printed (int precision) const { return size * (25+1); }
157 int Lvctr :: length_printed_vector (int precision) const { return 3 * (25+1); }
158 int Lvctr :: length_printed_tensor (int precision) const { return 9 * (25+1); }
159 
160 //* *** VECTOR
161 #ifdef DEBUG
162 void Lvctr :: printYourself (bool force) const
164 {
165  if (size > 20 && !force) _errorr("Count of rows is > 20, to force printing set first parameter as true");
166 
167  fprintf (stdout, "Long vector [%ld]:\n", size);
168 
169  for (long i=0; i<size; i++) fprintf (stdout, "%ld ", a[i]);
170 
171  fprintf (stdout, "\n");
172 }
173 #endif
174 
176 void Lvctr :: print (char *stream, int precision, double absZero) const
177 {
178  if (size < 1) errol;
179 
180  int numlen = precision+1;
181 
182  ; sprintf (stream, "%ld", a[0]); stream += numlen-1;
183  for (long i=1; i<size; i++) { sprintf (stream, " %*ld", precision, a[i]); stream += numlen; }
184 }
186 void Lvctr :: print_vector (char *stream, int precision, double absZero, bool zerorest) const
187 {
188  if (zerorest) { if (size < 0) errol; }
189  else if (size < 3) errol;
190 
191  sprintf (stream, "%*ld %*ld %*ld",
192  precision, size>=1 ? a[0] : 0,
193  precision, size>=2 ? a[1] : 0,
194  precision, size>=3 ? a[2] : 0);
195 }
196 
197 //* *** TENSOR
199 void pprint_symtensor (char *stream, int precision, const long *a)
200 {
201  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[0], precision, a[5], precision, a[4]); while (stream[0] != '\0') stream++;
202  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[5], precision, a[1], precision, a[3]); while (stream[0] != '\0') stream++;
203  sprintf (stream, "%*ld %*ld %*ld\n", precision, a[4], precision, a[3], precision, a[2]);
204 }
205 void pprint_symtensor_zeroL (char *stream, int precision)
206 {
207  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0); while (stream[0] != '\0') stream++;
208  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0); while (stream[0] != '\0') stream++;
209  sprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
210 }
211 void pprint_symtensor (FILE *stream, int precision, const long *a)
212 {
213  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[0], precision, a[5], precision, a[4]);
214  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[5], precision, a[1], precision, a[3]);
215  fprintf (stream, "%*ld %*ld %*ld\n", precision, a[4], precision, a[3], precision, a[2]);
216 }
217 void pprint_symtensor_zeroL (FILE *stream, int precision)
218 {
219  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
220  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
221  fprintf (stream, "%*d %*d %*d\n", precision, 0, precision, 0, precision, 0);
222 }
223 
224 void Lvctr :: print_symtensor (char *stream, int precision, double absZero, bool zerorest) const
225 {
226  if (zerorest) { if (size != 0 && size != 6) errol; }
227  else if (size != 6) errol;
228 
229  if (size) pprint_symtensor (stream, precision, a);
230  else pprint_symtensor_zeroL (stream, precision);
231 }
232 void Lvctr :: print_symtensor (FILE *stream, int precision, double absZero, bool zerorest) const
233 {
234  if (zerorest) { if (size != 0 && size != 6) errol; }
235  else if (size != 6) errol;
236 
237  if (size) pprint_symtensor (stream, precision, a);
238  else pprint_symtensor_zeroL (stream, precision);
239 }
240 
241 
242 
243 //* ********************************************************************************************
244 //* *** *** *** *** CLASS DVCTR *** *** *** ***
245 //* ********************************************************************************************
246 //* *** SET ***
248 void Dvctr :: realloc (long newsize)
249 {
250  if (newsize <= asize) return;
251 
252  asize = newsize;
253  double* newa = new double[asize];
254 
255  for (newsize=0; newsize<size; newsize++) *newa++ = *a++;
256 
257  delete [] (a-size);
258  a = newa-size;
259 }
260 
262 Dvctr* Dvctr :: assign_array (double *array, long s)
263 {
264  if (a) errol;
265  cnst = false;
266  exta = true;
267  size = s;
268  asize = 0;
269  a = array;
270  return this;
271 }
272 Dvctr* Dvctr :: assign_array (const double *array, long s)
273 {
274  if (a)
275  errol;
276  cnst = true;
277  exta = true;
278  size = s;
279  asize = 0;
280  a = (double *)array;
281  return this;
282 }
283 Dvctr* Dvctr :: free (void) {
284  if (a && !exta) delete [] a;
285  a = NULL;
286  exta = cnst = false;
287  shft = size = asize = 0;
288  return this;
289 }
290 
291 //* *** COMPUTE ***
292 double Dvctr :: give_lenght (void) const
293 {
294 # ifdef DEBUG
295  if (size == 0) _warningg ("size = 0 in give_length function");
296 # endif
297 
298  double length = 0.0;
299  for (long i=0; i<size; i++)
300  length += a[i] * a[i];
301 
302  return sqrt(length);
303 }
304 double Dvctr :: give_sum (void) const
305 {
306  double sum = 0;
307  for (long i=0; i<size; i++)
308  sum += a[i];
309 
310  return sum;
311 }
312 
313 //* *** EDIT ***
315 void Dvctr :: be_tnsr (const Dmtrx &src)
316 {
317  AVOID_CNST;
318  if ( src.give_crows() != 3 && src.give_ccols() != 3 ) _errorr ("size of mtrx != 3x3");
319 
320  const double *asrc = src.give_ptr2val(0,0);
321  this->resize_ignore_vals(6);
322 
323  a[0] = asrc[0];
324  a[1] = asrc[4];
325  a[2] = asrc[8];
326 
327  a[3] = asrc[5];
328  a[4] = asrc[2];
329  a[5] = asrc[1];
330 }
331 
333 void Dvctr :: be_vectproduct (const Dvctr &v1, const Dvctr &v2)
334 {
335  AVOID_CNST;
336  if ( v1.give_size() != 3 || v2.give_size() != 3 )
337  _errorr ("size of vectors is not equal to 3");
338 
339  this->resize_ignore_vals(3);
340 
341  const double *a1 = v1.give_ptr2val();
342  const double *a2 = v2.give_ptr2val();
343 
344  a[0] = a1[1] * a2[2] - a2[1] * a1[2];
345  a[1] = a1[2] * a2[0] - a2[2] * a1[0];
346  a[2] = a1[0] * a2[1] - a2[0] * a1[1];
347 }
348 
350 void Dvctr :: be_mean_of (const Dmtrx *src)
351 {
352  long n = src->give_crows();
353  long m = src->give_ccols();
354  const double *psrc = src->give_ptr2val();
355 
356  this->beCopyOf(psrc, m);
357 
358  if (n>1) {
359  psrc += m;
360  while (--n)
361  for (long i=0; i<m; i++)
362  a[i] += *(psrc++);
363 
364  this->dvdby(src->give_crows());
365  }
366 }
367 
369 double Dvctr :: give_dotproduct (const Dvctr &v)
370 {
371  double prd = 0.0;
372  for (long i=0; i<this->size; i++)
373  prd += a[i] * v.a[i];
374 
375  return prd;
376 }
377 
378 
380 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]; }
381 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]; }
382 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]; }
384 void Dvctr :: dvdby (double val)
385 {
386  AVOID_CNST;
387  for (long i=0; i<size; i++)
388  a[i] = a[i] / val;
389 }
391 void Dvctr :: tmsby (double val)
392 {
393  AVOID_CNST;
394  for (long i=0; i<size; i++)
395  a[i] = a[i] * val;
396 }
398 //void Dvctr :: tnsrRotWith1 (const Dmtrx &rot)
399 //{
400 // if ( size != 6 ) _errorr ("size of vectors is not equal to 6");
401 //
402 // Dmtrx sss, rotT;
403 //
404 // sss.be_tnsr(*this);
405 // rotT.be_transposition(rot);
406 //
407 // sss.tmsLeftBy (rot);
408 // sss.tmsRightBy (rotT);
409 //
410 // this->be_tnsr(sss);
411 //}
413 void Dvctr :: tnsrRotWith (const Dmtrx &rot)
414 {
415  AVOID_CNST;
416  if ( size != 6 ) _errorr ("size of vectors is not equal to 6");
417 
418  double aa[6];
419 
420  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];
421  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];
422  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];
423 
424  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];
425  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];
426  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];
427 
428  this->beCopyOf(aa);
429 }
430 
432 //void Dvctr :: tnsrRotAxisZangle (double a)
433 //{
434 // Dmtrx rot(3,3);
435 // rot(0,0) = cos(a); rot(0,1) = -sin(a); rot(0,2) = 0;
436 // rot(1,0) = sin(a); rot(1,1) = cos(a); rot(1,2) = 0;
437 // rot(2,0) = 0; rot(2,1) = 0; rot(2,2) = 1;
438 //
439 // this->tnsrRotWith(rot);
440 //}
442 void Dvctr :: tnsrRotAxisZangle (double alpha)
443 {
444  AVOID_CNST;
445  if ( size != 6 ) _errorr ("size of vectors is not equal to 6");
446 
447  double c = cos(alpha);
448  double s = sin(alpha);
449 
450  double aa[6];
451 
452  aa[0] = c*c * a[0] + s*s * a[1] - 2*c*s * a[5] ;
453  aa[1] = s*s * a[0] + c*c * a[1] + 2*s*c * a[5] ;
454  aa[2] = a[2] ;
455 
456  aa[5] = c*s * a[0] - s*c * a[1] + (c*c - s*s) * a[5] ;
457  aa[3] = c * a[3] + s * a[4];
458  aa[4] = -s * a[3] + c * a[4];
459 
460  this->beCopyOf(aa);
461 }
462 
463 
464 
466 void Dvctr :: normalize (void)
467 {
468  AVOID_CNST;
469  double length = this->give_lenght();
470  this->dvdby(length);
471 }
472 
474 //void Dvctr :: round2zero (double relZero, double absZero)
475 //{
476 // AVOID_CNST;
477 // double l;
478 // l = this->give_lenght();
479 // for (long i=0; i<size; i++)
480 // if (fabs(a[i]/l) < relZero || fabs(a[i]) < absZero) a[i] = 0.0;
481 //
482 //}
483 
484 
486 //void Dvctr :: round2zero (double relZero, double absZero)
487 //{
488 // AVOID_CNST;
489 // double l;
490 // l = this->give_lenght();
491 // for (long i=0; i<size; i++)
492 // if (fabs(a[i]/l) < relZero || fabs(a[i]) < absZero) a[i] = 0.0;
493 //
494 //}
495 
496 //if (zero < 0.0) {
497 // double l = 0.0;
498 // for (long i=0; i<n; i++) l += a[i] * a[i];
499 // l = sqrt(l) * zero;
500 //
501 // if ()
502 //
503 
504 
505 
506 
507 //* ***
508 //* *** PRINT ***
509 //* ***
510 int Dvctr :: length_printed (int precision) const { return size * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }
511 int Dvctr :: length_printed_vector (int precision) const { return 3 * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }
512 int Dvctr :: length_printed_tensor (int precision) const { return 9 * (1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT) + 1; }
513 
514 //* *** VECTOR
515 #ifdef DEBUG
516 void Dvctr :: printYourself (bool force) const
518 {
519  if (size > 20 && !force) _errorr("Count of rows is > 20, to force printing set first parameter as true");
520 
521  fprintf (stdout, "Double vector [%ld]:\n", size);
522 
523  for (long i=0; i<size; i++) fprintf (stdout, "% .*e ", PYP, a[i]);
524 
525  fprintf (stdout, "\n");
526 }
527 #endif
528 
530 void Dvctr :: print (char *stream, int precision, double absZero) const
531 {
532  if (size < 1) errol;
533 
534  double *p = a;
535  if (absZero != 0.0) p = give_copy_r2z (a, size, absZero);
536 
537  int numlen = 1+2+1+precision+2+NUM_DIGITS_IN_PRINTED_EXPONENT;
538 
539  ; sprintf (stream, "% .*e", precision, p[0]); stream += numlen-1;
540  for (long i=1; i<size; i++) { sprintf (stream, " % .*e", precision, p[i]); stream += numlen; }
541 
542  if (p != a) delete [] p;
543 }
545 void Dvctr :: print_vector (char *stream, int precision, double absZero, bool zerorest) const
546 {
547  if (zerorest) { if (size < 0) errol; }
548  else if (size < 3) errol;
549 
550  double *p = a;
551  if (absZero != 0.0 && size) p = give_copy_r2z (a, 3, absZero);
552 
553  sprintf (stream, "% .*e % .*e % .*e",
554  precision, size>=1 ? p[0] : 0.0,
555  precision, size>=2 ? p[1] : 0.0,
556  precision, size>=3 ? p[2] : 0.0);
557 
558  if (p != a) delete [] p;
559 }
560 
561 //* *** TENSOR
562 
564 void pprint_tensor (FILE *stream, int precision, const double *a)
565 {
566  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[1], precision, a[2]);
567  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[3], precision, a[4], precision, a[5]);
568  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[6], precision, a[7], precision, a[8]);
569 }
570 void pprint_symtensor (char *stream, int precision, const double *a)
571 {
572  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]); while (stream[0] != '\0') stream++;
573  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]); while (stream[0] != '\0') stream++;
574  sprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
575 }
576 void pprint_symtensor (FILE *stream, int precision, const double *a)
577 {
578  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);
579  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);
580  fprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
581 }
582 
583 void pprint_tensor_zero (char *stream, int precision)
584 {
585  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0); while (stream[0] != '\0') stream++;
586  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0); while (stream[0] != '\0') stream++;
587  sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
588 }
589 void pprint_tensor_zero (FILE *stream, int precision)
590 {
591  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
592  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
593  fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
594 }
595 
596 
597 
598 void Dvctr :: print_symtensor (char *stream, int precision, double absZero, bool zerorest) const
599 {
600  if (zerorest) { if (size != 0 && size != 6) errol; }
601  else if (size != 6) errol;
602 
603  double *p = a;
604  if (absZero != 0.0 && size) p = give_copy_r2z (a, 6, absZero);
605 
606  if (size) pprint_symtensor (stream, precision, p);
607  else pprint_tensor_zero (stream, precision);
608 
609  if (p != a) delete [] p;
610 }
611 void Dvctr :: print_symtensor (FILE *stream, int precision, double absZero, bool zerorest) const
612 {
613  if (zerorest) { if (size != 0 && size != 6) errol; }
614  else if (size != 6) errol;
615 
616  double *p = a;
617  if (absZero != 0.0 && size) p = give_copy_r2z (a, 6, absZero);
618 
619  if (size) pprint_symtensor (stream, precision, p);
620  else pprint_tensor_zero (stream, precision);
621 
622  if (p != a) delete [] p;
623 }
624 
625 
627 //void Dvctr :: print_symtensor (char *stream, int precision) const
628 //{
629 // sprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]); while (stream[0] != '\0') stream++;
630 // sprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]); while (stream[0] != '\0') stream++;
631 // sprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
632 //}
633 //void Dvctr :: print_symtensor_zero (char *stream, int precision) const
634 //{
635 // sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0); while (stream[0] != '\0') stream++;
636 // sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0); while (stream[0] != '\0') stream++;
637 // sprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
638 //}
640 //void Dvctr :: print_symtensor (FILE *stream, int precision) const
641 //{
642 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[5], precision, a[4]);
643 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[5], precision, a[1], precision, a[3]);
644 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[4], precision, a[3], precision, a[2]);
645 //}
646 //void Dvctr :: print_symtensor_zero (FILE *stream, int precision) const
647 //{
648 // fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
649 // fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
650 // fprintf (stream, "% .*e % .*e % .*e\n", precision, 0.0, precision, 0.0, precision, 0.0);
651 //}
652 
653 
654 
655 
656 
657 //* ********************************************************************************************
658 //* *** *** *** *** CLASS LMTRX *** *** *** ***
659 //* ********************************************************************************************
660 
661 //* *************
662 //* *** SET ***
663 //* *************
665 void Lmtrx :: resize_ignore_vals (long r, long c)
666 {
667  if (crow == r && ccol == c) return;
668 
669  crow = r; ccol = c;
670 
671  if (asize < crow*ccol) {
672  if (a) delete [] a;
673  asize = crow*ccol;
674  a = new long[asize];
675  }
676 }
677 
678 
679 //* ***************
680 //* *** PRINT ***
681 //* ***************
683 #ifdef DEBUG
684 void Lmtrx :: printYourself (bool force) const
686 {
687  if (crow > 50) _warningg("count of rows is > 50");
688  if (ccol > 50) _warningg("count of collumns is > 50");
689 
690  fprintf (stdout, "Long matrix [%ld,%ld]:\n", crow, ccol);
691 
692  long i, j;
693  for (i=0; i<crow; i++) {
694  for (j=0; j<ccol; j++)
695  fprintf (stdout, "%ld", a[i*ccol + j]);
696 
697  fprintf (stdout, "\n");
698  }
699 }
700 #endif
701 
702 
703 
704 //* ********************************************************************************************
705 //* *** *** *** *** CLASS DMTRX *** *** *** ***
706 //* ********************************************************************************************
707 
708 //* *************
709 //* *** SET ***
710 //* *************
711 void Dmtrx :: resize_ignore_vals (long r, long c)
712 {
713  if (crow == r && ccol == c) return;
714 
715  crow = r; ccol = c;
716 
717  if (asize < crow*ccol) {
718  if (a) delete [] a;
719  asize = crow*ccol;
720  a = new double[asize];
721  }
722 }
723 
724 
725 //* **************
726 //* *** EDIT ***
727 //* **************
729 void Dmtrx :: be_tnsr (const Dvctr &src)
730 {
731  if ( src.give_size() != 6 ) _errorr ("size of vectors is not equal to 6");
732 
733  const double *asrc = src.give_ptr2val();
734  this->resize_ignore_vals(3,3);
735 
736  a[0] = asrc[0];
737  a[4] = asrc[1];
738  a[8] = asrc[2];
739 
740  a[5] = a[7] = asrc[3];
741  a[2] = a[6] = asrc[4];
742  a[1] = a[3] = asrc[5];
743 }
745 void Dmtrx :: be_transposition (const Dmtrx &src)
746 {
747  this->resize_ignore_vals(src.give_ccols(), src.give_crows());
748 
749  const double *asrc = src.give_ptr2val(0,0);
750  long i, j;
751  for (i=0; i<crow; i++)
752  for (j=0; j<ccol; j++)
753  a[i*ccol + j] = asrc[j*ccol + i];
754 
755 }
757 void Dmtrx :: be_mtrxMmtrx (const Dmtrx &A, const Dmtrx &B)
758 {
759  if ( A.give_ccols() != B.give_crows() ) _errorr("size mismatch");
760  long ccr = A.give_ccols();
761 
762  this->resize_ignore_vals(A.give_crows(), B.give_ccols());
763  this->zero();
764 
765  const double *aA = A.give_ptr2val(0,0);
766  const double *aB = B.give_ptr2val(0,0);
767 
768  long i, j, k;
769  for (i=0; i<crow; i++)
770  for (j=0; j<ccol; j++)
771  for (k=0; k<ccr; k++)
772  a[i*ccol + j] += aA[i*ccr+k] * aB[k*ccol+j];
773 
774 }
776 void Dmtrx :: tmsLeftBy (const Dmtrx &src)
777 {
778  Dmtrx mtrxXthis;
779  mtrxXthis.be_mtrxMmtrx (src, *this);
780 
781  this->beCopyOf(mtrxXthis);
782 }
784 void Dmtrx :: tmsRightBy (const Dmtrx &src)
785 {
786  Dmtrx thisXmtrx;
787  thisXmtrx.be_mtrxMmtrx (*this, src);
788 
789  this->beCopyOf(thisXmtrx);
790 }
791 
792 
793 //* *****************
794 //* *** COMPUTE ***
795 //* *****************
797 double Dmtrx :: give_determinant (void) const
798 {
799  if (crow == ccol && ccol == 2) return a[0] * a[3] - a[1] * a[2];
800  else if (crow == ccol && ccol == 3) return a[0] * a[4] * a[8] + a[3] * a[7] * a[2] +
801  a[6] * a[1] * a[5] - a[6] * a[4] * a[2] -
802  a[7] * a[5] * a[0] - a[8] * a[3] * a[1] ;
803  else _errorr3("cannot compute determinant of matrix %ld x %ld", crow, ccol);
804  return 0.0;
805 }
806 
811 bool Dmtrx :: GaussSolve (Dvctr &x, const Dvctr &rhs)
812 {
813 #ifdef DEBUG
814  if (crow != ccol) _errorr3 ("matrix is not square, %ld != %ld", crow, ccol);
815  if (crow != rhs.give_size()) _errorr3 ("matrix and rhs dimension differs, %ld != %ld", crow, rhs.give_size());
816 #endif
817 
818  int i, j, k, maxrow;
819  double tmp;
820 
821  x = rhs;
822 
823  for (i=0; i<crow-1; i++) {
824 
825  // Find the row with the largest first value
826  maxrow = i;
827  for (j=i+1; j<crow; j++)
828  if (fabs(a[j*ccol+i]) > fabs(a[maxrow*ccol+i]))
829  maxrow = j;
830 
831  // Singular matrix?
832  if (fabs(a[maxrow*ccol+i]) < 1.0e-15)
833  return(false);
834 
835  // Swap the maxrow and ith row
836  if (maxrow != i) {
837  for (j=i; j<crow; j++) {
838  tmp = a[i*ccol+j];
839  a[i*ccol+j] = a[maxrow*ccol+j];
840  a[maxrow*ccol+j] = tmp;
841  }
842 
843  tmp = x[i];
844  x[i] = x[maxrow];
845  x[maxrow] = tmp;
846  }
847 
848  // Eliminate the ith element of the jth row
849  for (j=i+1; j<crow; j++) {
850  tmp = a[j*ccol+i] / a[i*ccol+i];
851 
852  for (k=i; k<crow; k++)
853  a[j*ccol+k] -= a[i*ccol+k] * tmp;
854 
855  x[j] -= x[i] * tmp;
856  }
857  }
858 
859  // Singular matrix?
860  if (fabs(a[crow*ccol-1]) < 1.0e-15)
861  return(false);
862 
863  // Do the back substitution
864  for (i=crow-1; i>=0; i--) {
865  tmp = 0;
866  for (j=i+1; j<crow; j++)
867  tmp += a[i*ccol+j] * x[j];
868 
869  x[i] = (x[i] - tmp) / a[i*ccol+i];
870  }
871 
872  return true;
873 }
874 
875 
880 bool Dmtrx :: GaussSolve (Dvctr *x, const Dvctr *rhs, long nrhs)
881 {
882 #ifdef DEBUG
883  if (crow != ccol) _errorr3 ("matrix is not square, %ld != %ld", crow, ccol);
884  if (crow != rhs[0].give_size()) _errorr3 ("matrix and rhs dimension differs, %ld != %ld", crow, rhs[0].give_size());
885 #endif
886 
887  int i, j, k, n, maxrow;
888  double tmp;
889 
890  for (n=0; n<nrhs; n++) x[n].beCopyOf(rhs[n]);
891 
892  for (i=0; i<crow-1; i++) {
893 
894  // Find the row with the largest first value
895  maxrow = i;
896  for (j=i+1; j<crow; j++)
897  if (fabs(a[j*ccol+i]) > fabs(a[maxrow*ccol+i]))
898  maxrow = j;
899 
900  // Singular matrix?
901  if (fabs(a[maxrow*ccol+i]) < 1.0e-15)
902  return(false);
903 
904  // Swap the maxrow and ith row
905  if (maxrow != i) {
906  for (j=i; j<crow; j++) {
907  tmp = a[i*ccol+j];
908  a[i*ccol+j] = a[maxrow*ccol+j];
909  a[maxrow*ccol+j] = tmp;
910  }
911 
912  for (n=0; n<nrhs; n++) {
913  tmp = x[n][i];
914  x[n][i] = x[n][maxrow];
915  x[n][maxrow] = tmp;
916  }
917  }
918 
919  // Eliminate the ith element of the jth row
920  for (j=i+1; j<crow; j++) {
921  tmp = a[j*ccol+i] / a[i*ccol+i];
922 
923  for (k=i; k<crow; k++)
924  a[j*ccol+k] -= a[i*ccol+k] * tmp;
925 
926  for (n=0; n<nrhs; n++)
927  x[n][j] -= x[n][i] * tmp;
928  }
929  }
930 
931  // Singular matrix?
932  if (fabs(a[crow*ccol-1]) < 1.0e-15)
933  return(false);
934 
935  // Do the back substitution
936  for (i=crow-1; i>=0; i--) {
937  for (n=0; n<nrhs; n++) {
938  tmp = 0;
939  for (j=i+1; j<crow; j++)
940  tmp += a[i*ccol+j] * x[n][j];
941 
942  x[n][i] = (x[n][i] - tmp) / a[i*ccol+i];
943  }
944  }
945 
946  return true;
947 }
948 
949 
950 //* ***************
951 //* *** PRINT ***
952 //* ***************
954 #ifdef DEBUG
955 void Dmtrx :: printYourself (bool force) const
957 {
958  if (crow > 10) _errorr("count of rows is > 10");
959  if (ccol > 10) _errorr("count of collumns is > 10");
960 
961  fprintf (stdout, "Double matrix [%ld,%ld]:\n", crow, ccol);
962 
963  long i, j;
964  for (i=0; i<crow; i++) {
965  for (j=0; j<ccol; j++)
966  fprintf (stdout, "% .*e ", PYP, a[i*ccol + j]);
967 
968  fprintf (stdout, "\n");
969  }
970 }
971 #endif
972 
973 
974 //void Dmtrx :: print_tensor (FILE *stream, int precision) const
975 //{
976 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[0], precision, a[1], precision, a[2]);
977 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[3], precision, a[4], precision, a[5]);
978 // fprintf (stream, "% .*e % .*e % .*e\n", precision, a[6], precision, a[7], precision, a[8]);
979 //}
980 
981 void Dmtrx :: print_tensor (FILE *stream, int precision, double absZero, bool zerorest) const
982 {
983  if (zerorest) { if (! ((crow == 0 && ccol == 0) || (crow == 3 && ccol == 3)) ) errol; }
984  else if (crow != 3 || ccol != 3) errol;
985 
986  double *p = a;
987  if (absZero != 0.0 && crow) p = give_copy_r2z (a, 9, absZero);
988 
989  if (crow) pprint_tensor (stream, precision, p);
990  else pprint_tensor_zero (stream, precision);
991 
992  if (p != a) delete [] p;
993 }
994 
995 
996 } // namespace midaspace
midaspace::Lvctr::free
Lvctr * free(void)
Definition: arrays.cpp:131
midaspace::Dmtrx::beCopyOf
void beCopyOf(const double *src)
Definition: arrays.h:808
_warningg
#define _warningg(_1)
Definition: gelib.h:154
midaspace::Xvctr::give_lenght
virtual double give_lenght(void) const =0
midaspace::Lvctr::a
long * a
Array of values.
Definition: arrays.h:411
midaspace::Lvctr::length_printed_vector
int length_printed_vector(int precision) const
Definition: arrays.cpp:157
midaspace::Dmtrx::GaussSolve
bool GaussSolve(Dvctr &x, const Dvctr &rhs)
gaussian eliminatio
Definition: arrays.cpp:811
midaspace::Dvctr::length_printed
int length_printed(int precision) const
void round2zero (double zero, double absZero);
Definition: arrays.cpp:510
midaspace::Dmtrx::tmsLeftBy
void tmsLeftBy(const Dmtrx &src)
Definition: arrays.cpp:776
midaspace::Xvctr::give_zero
double give_zero(double abszero, double relzero) const
Definition: arrays.cpp:85
midaspace::Lvctr::length_printed
int length_printed(int precision) const
Definition: arrays.cpp:156
gelib.h
General functions.
midaspace::Dvctr::dvdby
void dvdby(double val)
Definition: arrays.cpp:384
midaspace::Dvctr::give_dotproduct
double give_dotproduct(const Dvctr &v)
Definition: arrays.cpp:369
midaspace::Dvctr::length_printed_tensor
int length_printed_tensor(int precision) const
Definition: arrays.cpp:512
midaspace::Dmtrx::print_tensor
void print_tensor(FILE *stream, int precision, double absZero, bool zerorest) const
print yourself
Definition: arrays.cpp:981
midaspace::Dvctr::realloc
void realloc(long newsize)
reallocate up receiver
Definition: arrays.cpp:248
midaspace::Lvctr::print_vector
void print_vector(char *stream, int precision, double absZero=0.0, bool zerorest=false) const
Definition: arrays.cpp:186
midaspace::Dmtrx::be_mtrxMmtrx
void be_mtrxMmtrx(const Dmtrx &A, const Dmtrx &B)
Definition: arrays.cpp:757
midaspace::PoinT::beRotatedPoint
void beRotatedPoint(const Dmtrx *rot, const PoinT *point)
Receiver will be point 'point' rotated by matrix 'rot'. this = rot . point.
Definition: arrays.cpp:60
midaspace::Xvctr::give_size
virtual long give_size(void) const
Definition: arrays.h:387
midaspace::Dvctr::normalize
void normalize(void)
Definition: arrays.cpp:466
midaspace::Lvctr::realloc
void realloc(long newsize)
reallocate up receiver
Definition: arrays.cpp:97
NUM_DIGITS_IN_PRINTED_EXPONENT
#define NUM_DIGITS_IN_PRINTED_EXPONENT
Definition: gelib.h:31
midaspace::Dmtrx::tmsRightBy
void tmsRightBy(const Dmtrx &src)
Definition: arrays.cpp:784
midaspace::Lvctr::assign_array
Lvctr * assign_array(long *array, long s)
Definition: arrays.cpp:111
midaspace::pprint_tensor_zero
void pprint_tensor_zero(char *stream, int precision)
Definition: arrays.cpp:583
errol
#define errol
Definition: gelib.h:142
midaspace::pprint_symtensor_zeroL
void pprint_symtensor_zeroL(char *stream, int precision)
Definition: arrays.cpp:205
midaspace::give_copy_r2z
double give_copy_r2z(double a, double zero)
Definition: arrays.cpp:13
midaspace::Dvctr::give_sum
double give_sum(void) const
Definition: arrays.cpp:304
midaspace::Dmtrx::give_determinant
double give_determinant(void) const
Definition: arrays.cpp:797
midaspace::Dmtrx::be_transposition
void be_transposition(const Dmtrx &src)
Definition: arrays.cpp:745
midaspace::Dmtrx::resize_ignore_vals
void resize_ignore_vals(long r, long c)
print yourself
Definition: arrays.cpp:711
midaspace::Dvctr::sbt
void sbt(const Dvctr &src)
Definition: arrays.cpp:382
AVOID_CNST
#define AVOID_CNST
Definition: arrays.h:264
midaspace::Dvctr::be_vectproduct
void be_vectproduct(const Dvctr &v1, const Dvctr &v2)
Definition: arrays.cpp:333
arrays.h
Structs Elem3D, PoinT and VectoR; classes Array, Array1d, Xscal, Dscal, Xvctr, Lvctr, Dvctr, Xmtrx, Lmtrx and Dmtrx.
midaspace::Dvctr::resize_ignore_vals
Dvctr * resize_ignore_vals(long newsize)
Definition: arrays.h:595
midaspace::Dvctr::beCopyOf
void beCopyOf(const PoinT &src)
Definition: arrays.h:632
midaspace::Lvctr::give_sum
long give_sum(void) const
Definition: arrays.cpp:148
midaspace::Dvctr::print_vector
void print_vector(char *stream, int precision, double absZero=0.0, bool zerorest=false) const
Definition: arrays.cpp:545
midaspace::Dvctr::print
void print(char *stream, int precision, double absZero=0.0) const
Definition: arrays.cpp:530
midaspace::Dmtrx::be_tnsr
void be_tnsr(const Dvctr &src)
Definition: arrays.cpp:729
midaspace::Dvctr::give_lenght
double give_lenght(void) const
Definition: arrays.cpp:292
midaspace::Dvctr::free
Dvctr * free(void)
Definition: arrays.cpp:283
midaspace::Lvctr::print
void print(char *stream, int precision, double absZero=0.0) const
Definition: arrays.cpp:176
midaspace::Lmtrx::resize_ignore_vals
void resize_ignore_vals(long r, long c)
Definition: arrays.cpp:665
midaspace::Lvctr::give_lenght
double give_lenght(void) const
Definition: arrays.cpp:140
midaspace::Dvctr::add
void add(const Dvctr &src)
Definition: arrays.cpp:380
midaspace::MatriX::rotate
void rotate(Dvctr *v) const
Definition: arrays.cpp:30
_errorr
#define _errorr(_1)
Definition: gelib.h:151
midaspace::Dmtrx::zero
void zero(void)
Definition: arrays.h:799
_errorr3
#define _errorr3(_1, _2, _3)
Definition: gelib.h:146
midaspace::Dvctr::assign_array
Dvctr * assign_array(double *array, long s)
Definition: arrays.cpp:262
midaspace::Dscal::print
virtual void print(char *stream, int precision, double absZero=0.0) const
Definition: arrays.cpp:77
midaspace::Dvctr::be_mean_of
void be_mean_of(const Dmtrx *src)
Definition: arrays.cpp:350
midaspace::MatriX::a
double a[9]
Definition: arrays.h:226
midaspace::pprint_tensor
void pprint_tensor(FILE *stream, int precision, const double *a)
Definition: arrays.cpp:564
midaspace::Lvctr::length_printed_tensor
int length_printed_tensor(int precision) const
Definition: arrays.cpp:158
midaspace::Dvctr::a
double * a
Definition: arrays.h:559
midaspace::Lvctr::give_number_of_nonzeros
long give_number_of_nonzeros(void) const
Definition: arrays.cpp:149
midaspace::Dvctr::tnsrRotAxisZangle
void tnsrRotAxisZangle(double a)
rotate coord system around axis z by angle
Definition: arrays.cpp:442
midaspace::Dvctr::tnsrRotWith
void tnsrRotWith(const Dmtrx &rot)
this = rot * this * rot^transp
Definition: arrays.cpp:413
midaspace::Dvctr::be_tnsr
void be_tnsr(const Dmtrx &src)
Definition: arrays.cpp:315
midaspace::Dmtrx::give_ptr2val
double * give_ptr2val(long r=0, long c=0)
Definition: arrays.h:793
midaspace::Lvctr::give_number_of_zeros
long give_number_of_zeros(void) const
Definition: arrays.cpp:150
midaspace::Dmtrx::a
double * a
Definition: arrays.h:774
midaspace::Dvctr::give_ptr2val
double * give_ptr2val(long i=0)
return pointer to
Definition: arrays.h:618
midaspace::Dvctr::addtms
void addtms(const Dvctr &src, double tms)
Definition: arrays.cpp:381
midaspace::Dvctr::print_symtensor
void print_symtensor(char *stream, int precision, double absZero=0.0, bool zerorest=false) const
Definition: arrays.cpp:598
midaspace::Lvctr::print_symtensor
void print_symtensor(char *stream, int precision, double absZero=0.0, bool zerorest=false) const
Definition: arrays.cpp:224
midaspace::Xmtrx::give_ccols
long give_ccols(void) const
Definition: arrays.h:729
PYP
#define PYP
Definition: arrays.h:21
midaspace::Dvctr::length_printed_vector
int length_printed_vector(int precision) const
Definition: arrays.cpp:511
midaspace::Xmtrx::give_crows
long give_crows(void) const
Definition: arrays.h:728
midaspace::Dvctr::tmsby
void tmsby(double val)
Definition: arrays.cpp:391
midaspace::pprint_symtensor
void pprint_symtensor(char *stream, int precision, const long *a)
Definition: arrays.cpp:199
This reference manual is part of the MIDAS documentation. Copyright (c) 2005-2013 Ladislav Termit Svoboda. Generated on Fri Mar 21 2014 11:39:42 for MIDAS by doxygen 1.8.6.