SAGA API  v9.6
mat_matrix.cpp
Go to the documentation of this file.
1 
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // mat_matrix.cpp //
15 // //
16 // Copyright (C) 2005 by Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 // e-mail: oconrad@saga-gis.org //
49 // //
51 
52 //---------------------------------------------------------
53 #include "mat_tools.h"
54 
55 
57 // //
58 // //
59 // //
61 
62 //---------------------------------------------------------
65 
66 
68 // //
69 // //
70 // //
72 
73 //---------------------------------------------------------
75 {
76  m_Array.Create(sizeof(double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_2);
77 }
78 
79 //---------------------------------------------------------
81 {
82  m_Array.Create(sizeof(double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_2);
83 
84  Create(Vector);
85 }
86 
87 bool CSG_Vector::Create(const CSG_Vector &Vector)
88 {
89  return( Assign(Vector) );
90 }
91 
92 //---------------------------------------------------------
93 CSG_Vector::CSG_Vector(sLong n, const double *Data)
94 {
95  m_Array.Create(sizeof(double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_2);
96 
97  Create(n, Data);
98 }
99 
100 bool CSG_Vector::Create(sLong n, const double *Data)
101 {
102  if( n > 0 )
103  {
104  if( m_Array.Set_Array(n) )
105  {
106  if( Data )
107  {
108  memcpy(Get_Data(), Data, n * sizeof(double));
109  }
110  else
111  {
112  memset(Get_Data(), 0, n * sizeof(double));
113  }
114 
115  return( true );
116  }
117  }
118 
119  Destroy();
120 
121  return( false );
122 }
123 
124 //---------------------------------------------------------
126 {
127  Destroy();
128 }
129 
131 {
132  return( m_Array.Set_Array(0) );
133 }
134 
135 //---------------------------------------------------------
141 {
142  if( nRows > Get_Size() )
143  {
144  return( Add_Rows(nRows - Get_Size()) );
145  }
146 
147  if( nRows < Get_Size() )
148  {
149  return( Del_Rows(Get_Size() - nRows) );
150  }
151 
152  return( true );
153 }
154 
155 //---------------------------------------------------------
157 {
158  if( nRows > 0 && m_Array.Set_Array(Get_Size() + nRows) )
159  {
160  for(sLong i=Get_Size()-nRows; i<Get_Size(); i++)
161  {
162  Get_Data()[i] = 0.;
163  }
164 
165  return( true );
166  }
167 
168  return( false );
169 }
170 
171 //---------------------------------------------------------
178 {
179  if( nRows >= Get_Size() )
180  {
181  return( Destroy() );
182  }
183 
184  return( nRows < 1 ? false : m_Array.Set_Array(Get_Size() - nRows) );
185 }
186 
187 //---------------------------------------------------------
188 bool CSG_Vector::Add_Row(double Value)
189 {
190  if( m_Array.Inc_Array() )
191  {
192  Get_Data()[Get_Size() - 1] = Value;
193 
194  return( true );
195  }
196 
197  return( false );
198 }
199 
200 //---------------------------------------------------------
202 {
203  if( Row < Get_Size() )
204  {
205  if( Row >= 0 ) // just remove last entry => if( Row < 0 )
206  {
207  for(sLong i=Row, j=Row+1; j<Get_Size(); i++, j++)
208  {
209  Get_Data()[i] = Get_Data()[j];
210  }
211  }
212 
213  return( m_Array.Dec_Array() );
214  }
215 
216  return( false );
217 }
218 
219 
221 // //
223 
224 //---------------------------------------------------------
225 CSG_String CSG_Vector::to_String(int Width, int Precision, bool bScientific, const SG_Char *Separator) const
226 {
227  CSG_String s, sep(Separator && *Separator ? Separator : SG_T(" "));
228 
229  for(sLong i=0; i<Get_Size(); i++)
230  {
231  s += sep + SG_Get_Double_asString(Get_Data(i), Width, Precision, bScientific);
232  }
233 
234  return( s );
235 }
236 
237 //---------------------------------------------------------
239 {
240  Destroy();
241 
242  CSG_String_Tokenizer Line(String);
243 
244  while( Line.Has_More_Tokens() )
245  {
246  double d; CSG_String s(Line.Get_Next_Token());
247 
248  if( s.asDouble(d) )
249  {
250  Add_Row(d);
251  }
252  }
253 
254  return( Get_Size() > 0 );
255 }
256 
257 
259 // //
261 
262 //---------------------------------------------------------
263 bool CSG_Vector::is_Null(void) const
264 {
265  for(sLong i=0; i<Get_N(); i++)
266  {
267  if( Get_Data(i) != 0. )
268  {
269  return( false );
270  }
271  }
272 
273  return( true );
274 }
275 
276 //---------------------------------------------------------
277 bool CSG_Vector::is_Equal(const CSG_Vector &Vector) const
278 {
279  if( Get_Size() == Vector.Get_Size() )
280  {
281  for(sLong i=0; i<Get_N(); i++)
282  {
283  if( Get_Data(i) != Vector.Get_Data(i) )
284  {
285  return( false );
286  }
287  }
288 
289  return( true );
290  }
291 
292  return( false );
293 }
294 
295 //---------------------------------------------------------
296 bool CSG_Vector::is_Collinear(const CSG_Vector &Vector) const
297 {
298  if( Get_Size() == Vector.Get_Size() && !is_Null() && !Vector.is_Null() )
299  {
300  double b = Vector.Get_Length() / Get_Length();
301 
302  for(sLong i=0; i<Get_N(); i++)
303  {
304  if( b * Get_Data(i) != Vector.Get_Data(i) )
305  {
306  return( false );
307  }
308  }
309 
310  return( true );
311  }
312 
313  return( false );
314 }
315 
316 //---------------------------------------------------------
317 bool CSG_Vector::Assign(double Scalar)
318 {
319  if( Get_Size() > 0 )
320  {
321  for(sLong i=0; i<Get_Size(); i++)
322  {
323  Get_Data()[i] = Scalar;
324  }
325 
326  return( true );
327  }
328 
329  return( false );
330 }
331 
332 //---------------------------------------------------------
333 bool CSG_Vector::Assign(const CSG_Vector &Vector)
334 {
335  if( Create(Vector.Get_Size()) )
336  {
337  memcpy(Get_Data(), Vector.Get_Data(), Get_Size() * sizeof(double));
338 
339  return( true );
340  }
341 
342  return( false );
343 }
344 
345 //---------------------------------------------------------
346 bool CSG_Vector::Add(double Scalar)
347 {
348  if( Get_Size() > 0 )
349  {
350  for(sLong i=0; i<Get_Size(); i++)
351  {
352  Get_Data()[i] += Scalar;
353  }
354 
355  return( true );
356  }
357 
358  return( false );
359 }
360 
361 //---------------------------------------------------------
362 bool CSG_Vector::Add(const CSG_Vector &Vector)
363 {
364  if( Get_Size() == Vector.Get_Size() && Get_Size() > 0 )
365  {
366  for(sLong i=0; i<Get_Size(); i++)
367  {
368  Get_Data()[i] += Vector.Get_Data()[i];
369  }
370 
371  return( true );
372  }
373 
374  return( false );
375 }
376 
377 //---------------------------------------------------------
378 bool CSG_Vector::Subtract(const CSG_Vector &Vector)
379 {
380  if( Get_Size() == Vector.Get_Size() && Get_Size() > 0 )
381  {
382  for(sLong i=0; i<Get_Size(); i++)
383  {
384  Get_Data()[i] -= Vector.Get_Data()[i];
385  }
386 
387  return( true );
388  }
389 
390  return( false );
391 }
392 
393 //---------------------------------------------------------
394 bool CSG_Vector::Multiply(double Scalar)
395 {
396  if( Get_Size() > 0 )
397  {
398  for(sLong i=0; i<Get_Size(); i++)
399  {
400  Get_Data()[i] *= Scalar;
401  }
402 
403  return( true );
404  }
405 
406  return( false );
407 }
408 
409 //---------------------------------------------------------
410 bool CSG_Vector::Multiply(const CSG_Vector &Vector)
411 {
412  return( Multiply_Cross(Vector) );
413 }
414 
415 //---------------------------------------------------------
417 {
418  if( Get_Size() == Vector.Get_Size() && Get_Size() == 3 )
419  {
420  CSG_Vector v(*this);
421 
422  Get_Data()[0] = v[1] * Vector[2] - v[2] * Vector[1];
423  Get_Data()[1] = v[2] * Vector[0] - v[0] * Vector[2];
424  Get_Data()[2] = v[0] * Vector[1] - v[1] * Vector[0];
425 
426  return( true );
427  }
428 
429  return( false );
430 }
431 
432 //---------------------------------------------------------
433 double CSG_Vector::Multiply_Scalar(const CSG_Vector &Vector) const
434 {
435  double z = 0.;
436 
437  if( Get_Size() == Vector.Get_Size() )
438  {
439  for(sLong i=0; i<Get_Size(); i++)
440  {
441  z += Get_Data()[i] * Vector[i];
442  }
443  }
444 
445  return( z );
446 }
447 
448 //---------------------------------------------------------
449 bool CSG_Vector::Multiply(const class CSG_Matrix &Matrix)
450 {
451  return( Assign(Matrix.Multiply(*this)) );
452 }
453 
454 
456 // //
458 
459 //---------------------------------------------------------
461 {
462  Assign(Scalar);
463 
464  return( *this );
465 }
466 
468 {
469  Assign(Vector);
470 
471  return( *this );
472 }
473 
474 //---------------------------------------------------------
476 {
477  Add(Scalar);
478 
479  return( *this );
480 }
481 
483 {
484  Add(Vector);
485 
486  return( *this );
487 }
488 
489 //---------------------------------------------------------
491 {
492  Add(-Scalar);
493 
494  return( *this );
495 }
496 
498 {
499  Subtract(Vector);
500 
501  return( *this );
502 }
503 
504 //---------------------------------------------------------
506 {
507  Multiply(Scalar);
508 
509  return( *this );
510 }
511 
513 {
514  Multiply(Vector);
515 
516  return( *this );
517 }
518 
520 {
521  Multiply(Matrix);
522 
523  return( *this );
524 }
525 
526 //---------------------------------------------------------
528 {
529  CSG_Vector v(*this);
530 
531  v.Add(Scalar);
532 
533  return( v );
534 }
535 
537 {
538  CSG_Vector v(*this);
539 
540  v.Add(Vector);
541 
542  return( v );
543 }
544 
545 //---------------------------------------------------------
547 {
548  CSG_Vector v(*this);
549 
550  v.Add(-Scalar);
551 
552  return( v );
553 }
554 
556 {
557  CSG_Vector v(*this);
558 
559  v.Subtract(Vector);
560 
561  return( v );
562 }
563 
564 //---------------------------------------------------------
566 {
567  CSG_Vector v(*this);
568 
569  v.Multiply(Scalar);
570 
571  return( v );
572 }
573 
574 double CSG_Vector::operator * (const CSG_Vector &Vector) const
575 {
576  return( Multiply_Scalar(Vector) );
577 }
578 
579 CSG_Vector operator * (double Scalar, const CSG_Vector &Vector)
580 {
581  return( Vector * Scalar );
582 }
583 
584 
586 // //
588 
589 //---------------------------------------------------------
591 {
592  return( Create(Get_Size()) );
593 }
594 
595 //---------------------------------------------------------
597 {
598  double Length = Get_Length();
599 
600  if( Length > 0. )
601  {
602  for(sLong i=0; i<Get_Size(); i++)
603  {
604  Get_Data()[i] /= Length;
605  }
606 
607  return( true );
608  }
609 
610  return( false );
611 }
612 
613 
615 // //
617 
618 //---------------------------------------------------------
620 {
621  if( Get_Size() > 1 )
622  {
623  double *v = Get_Data();
624 
625  for(sLong i=0, j=Get_Size()-1; i<j; i++, j--)
626  {
627  double d = v[i]; v[i] = v[j]; v[j] = d;
628  }
629 
630  return( true );
631  }
632 
633  return( false );
634 }
635 
636 
638 // //
640 
641 //---------------------------------------------------------
642 bool CSG_Vector::Sort(bool bAscending)
643 {
644  if( Get_Size() > 0 )
645  {
646  qsort(Get_Data(), Get_Size(), sizeof(double), SG_Compare_Double);
647 
648  if( bAscending == false )
649  {
650  Flip_Values();
651  }
652 
653  return( true );
654  }
655 
656  return( false );
657 }
658 
659 
661 // //
663 
664 //---------------------------------------------------------
666 {
667  CSG_Vector v(*this);
668 
669  v.Multiply_Scalar(Vector);
670 
671  return( v );
672 }
673 
674 //---------------------------------------------------------
676 {
677  CSG_Vector v(*this);
678 
679  v.Multiply_Cross(Vector);
680 
681  return( v );
682 }
683 
684 //---------------------------------------------------------
685 double CSG_Vector::Get_Length(void) const
686 {
687  if( Get_Size() > 0 )
688  {
689  double z = 0., *Z = Get_Data();
690 
691  for(sLong i=0; i<Get_Size(); i++)
692  {
693  z += Z[i] * Z[i];
694  }
695 
696  return( sqrt(z) );
697  }
698 
699  return( 0. );
700 }
701 
702 //---------------------------------------------------------
703 double CSG_Vector::Get_Angle(const CSG_Vector &Vector) const
704 {
705  if( Get_Size() > Vector.Get_Size() )
706  {
707  return( Vector.Get_Angle(*this) );
708  }
709 
710  double A, B;
711 
712  if( (A = Get_Length()) > 0. && (B = Vector.Get_Length()) > 0. )
713  {
714  double z = 0., *Z = Get_Data();
715 
716  for(sLong i=0; i<Get_Size(); i++)
717  {
718  z += Vector[i] * Z[i];
719  }
720 
721  for(sLong i=Get_Size(); i<Vector.Get_Size(); i++)
722  {
723  z += Vector[i];
724  }
725 
726  return( acos(z / (A * B)) );
727  }
728 
729  return( 0. );
730 }
731 
732 //---------------------------------------------------------
734 {
735  CSG_Vector v(*this);
736 
737  v.Set_Unity();
738 
739  return( v );
740 }
741 
742 
744 // //
745 // //
746 // //
748 
749 //---------------------------------------------------------
750 bool SG_VectorR2_Rotate(double &x, double &y, double Angle)
751 {
752  double sin_a = sin(Angle);
753  double cos_a = cos(Angle);
754 
755  double t = x;
756 
757  x = t * cos_a - y * sin_a;
758  y = t * sin_a + y * cos_a;
759 
760  return( true );
761 }
762 
763 //---------------------------------------------------------
764 bool SG_VectorR2_Rotate(double Vector[2], double Angle)
765 {
766  return( SG_VectorR2_Rotate(Vector[0], Vector[1], Angle) );
767 }
768 
769 //---------------------------------------------------------
770 bool SG_VectorR2_Rotate(CSG_Vector &Vector , double Angle)
771 {
772  return( Vector.Get_Size() >= 2 && SG_VectorR2_Rotate(Vector[0], Vector[1], Angle) );
773 }
774 
775 //---------------------------------------------------------
776 bool SG_VectorR3_Rotate(double Vector[3], size_t Axis, double Angle)
777 {
778  if( Axis <= 3 )
779  {
780  CSG_Vector v(3, Vector);
781 
782  double sin_a = sin(Angle);
783  double cos_a = cos(Angle);
784 
785  switch( Axis )
786  {
787  case 0:
788  Vector[1] = v[1] * cos_a - v[2] * sin_a;
789  Vector[2] = v[1] * sin_a + v[2] * cos_a;
790  break;
791 
792  case 1:
793  Vector[0] = v[0] * cos_a + v[2] * sin_a;
794  Vector[2] =-v[0] * sin_a + v[2] * cos_a;
795  break;
796 
797  case 2:
798  Vector[0] = v[0] * cos_a - v[1] * sin_a;
799  Vector[1] = v[0] * sin_a + v[1] * cos_a;
800  break;
801  }
802 
803  return( true );
804  }
805 
806  return( false );
807 }
808 
809 //---------------------------------------------------------
810 bool SG_VectorR3_Rotate(CSG_Vector &Vector , size_t Axis, double Angle)
811 {
812  return( Vector.Get_Size() >= 3 && SG_VectorR3_Rotate(Vector.Get_Data(), Axis, Angle) );
813 }
814 
815 
817 // //
818 // //
819 // //
821 
822 //---------------------------------------------------------
824 {
825  _On_Construction();
826 }
827 
828 //---------------------------------------------------------
830 {
831  _On_Construction();
832 
833  Create(Matrix);
834 }
835 
836 bool CSG_Matrix::Create(const CSG_Matrix &Matrix)
837 {
838  return( Create(Matrix.m_nx, Matrix.m_ny, Matrix.m_z[0]) );
839 }
840 
841 //---------------------------------------------------------
842 CSG_Matrix::CSG_Matrix(sLong nCols, sLong nRows, const double *Data)
843 {
844  _On_Construction();
845 
846  Create(nCols, nRows, Data);
847 }
848 
849 bool CSG_Matrix::Create(sLong nCols, sLong nRows, const double *Data)
850 {
851  if( nCols < 1 || nRows < 1 )
852  {
853  Destroy();
854 
855  return( false );
856  }
857 
858  if( nCols != m_nx || nRows != m_ny )
859  {
860  Destroy();
861 
862  if( (m_z = (double **)SG_Malloc(nRows * sizeof(double *))) == NULL
863  || (m_z[0] = (double *)SG_Malloc(nRows * nCols * sizeof(double ))) == NULL )
864  {
865  Destroy();
866 
867  return( false );
868  }
869 
870  m_nx = nCols;
871  m_ny = nRows;
872 
873  for(sLong y=1; y<m_ny; y++)
874  {
875  m_z[y] = m_z[y - 1] + m_nx;
876  }
877  }
878 
879  if( Data )
880  {
881  memcpy(m_z[0], Data, m_ny * m_nx * sizeof(double));
882  }
883  else
884  {
885  memset(m_z[0], 0, m_ny * m_nx * sizeof(double));
886  }
887 
888  return( true );
889 }
890 
891 //---------------------------------------------------------
892 CSG_Matrix::CSG_Matrix(sLong nCols, sLong nRows, const double **Data)
893 {
894  _On_Construction();
895 
896  Create(nCols, nRows, Data);
897 }
898 
899 bool CSG_Matrix::Create(sLong nCols, sLong nRows, const double **Data)
900 {
901  if( Create(nCols, nRows) )
902  {
903  if( Data )
904  {
905  for(sLong y=0; y<m_ny; y++)
906  {
907  memcpy(m_z[y], Data[y], m_nx * sizeof(double));
908  }
909  }
910 
911  return( true );
912  }
913 
914  return( false );
915 }
916 
917 //---------------------------------------------------------
919 {
920  Destroy();
921 }
922 
924 {
925  if( m_z )
926  {
927  if( m_z[0] )
928  {
929  SG_Free(m_z[0]);
930  }
931 
932  SG_Free(m_z);
933  }
934 
935  _On_Construction();
936 
937  return( true );
938 }
939 
940 //---------------------------------------------------------
941 void CSG_Matrix::_On_Construction(void)
942 {
943  m_z = NULL; m_nx = m_ny = 0;
944 }
945 
946 
948 // //
950 
951 //---------------------------------------------------------
952 bool CSG_Matrix::Set_Size(sLong nRows, sLong nCols)
953 {
954  return( nRows > 0 && nCols > 0 && Set_Rows(nRows) && Set_Cols(nCols) );
955 }
956 
957 //---------------------------------------------------------
959 {
960  if( nCols > m_nx )
961  {
962  return( Add_Cols(nCols - m_nx) );
963  }
964 
965  if( nCols < m_nx )
966  {
967  return( Del_Cols(m_nx - nCols) );
968  }
969 
970  return( true );
971 }
972 
973 //---------------------------------------------------------
975 {
976  if( nRows > m_ny )
977  {
978  return( Add_Rows(nRows - m_ny) );
979  }
980 
981  if( nRows < m_ny )
982  {
983  return( Del_Rows(m_ny - nRows) );
984  }
985 
986  return( true );
987 }
988 
989 //---------------------------------------------------------
991 {
992  if( nCols > 0 && m_ny > 0 )
993  {
994  CSG_Matrix Tmp(*this);
995 
996  if( Create(Tmp.m_nx + nCols, Tmp.m_ny) )
997  {
998  for(sLong y=0; y<Tmp.m_ny; y++)
999  {
1000  memcpy(m_z[y], Tmp.m_z[y], Tmp.m_nx * sizeof(double));
1001  }
1002 
1003  return( true );
1004  }
1005  }
1006 
1007  return( false );
1008 }
1009 
1010 //---------------------------------------------------------
1012 {
1013  if( nRows > 0 && m_nx > 0 )
1014  {
1015  m_ny += nRows;
1016 
1017  m_z = (double **)SG_Realloc(m_z , m_ny * sizeof(double *));
1018  m_z[0] = (double *)SG_Realloc(m_z[0], m_ny * m_nx * sizeof(double ));
1019 
1020  for(sLong y=1; y<m_ny; y++)
1021  {
1022  m_z[y] = m_z[y - 1] + m_nx;
1023  }
1024 
1025  memset(m_z[m_ny - nRows], 0, nRows * m_nx * sizeof(double));
1026 
1027  return( true );
1028  }
1029 
1030  return( false );
1031 }
1032 
1033 //---------------------------------------------------------
1038 {
1039  if( nCols > 0 && m_ny > 0 && nCols < m_nx )
1040  {
1041  CSG_Matrix Tmp(*this);
1042 
1043  if( Create(Tmp.m_nx - nCols, Tmp.m_ny) )
1044  {
1045  for(sLong y=0; y<Tmp.m_ny; y++)
1046  {
1047  memcpy(m_z[y], Tmp.m_z[y], m_nx * sizeof(double));
1048  }
1049 
1050  return( true );
1051  }
1052  }
1053 
1054  return( false );
1055 }
1056 
1057 //---------------------------------------------------------
1062 {
1063  if( nRows > 0 && m_nx > 0 && nRows < m_ny )
1064  {
1065  m_ny -= nRows;
1066 
1067  m_z = (double **)SG_Realloc(m_z , m_ny * sizeof(double *));
1068  m_z[0] = (double *)SG_Realloc(m_z[0], m_ny * m_nx * sizeof(double ));
1069 
1070  return( true );
1071  }
1072 
1073  return( false );
1074 }
1075 
1076 //---------------------------------------------------------
1077 bool CSG_Matrix::Add_Col(const double *Data)
1078 {
1079  if( Add_Cols(1) )
1080  {
1081  Set_Col(m_nx - 1, Data);
1082 
1083  return( true );
1084  }
1085 
1086  return( false );
1087 }
1088 
1090 {
1091  if( m_nx == 0 )
1092  {
1093  return( Create(1, Data.Get_Size()) && Set_Col(0, Data.Get_Data()) );
1094  }
1095 
1096  return( m_ny <= Data.Get_Size() && Add_Col(Data.Get_Data()) );
1097 }
1098 
1099 //---------------------------------------------------------
1100 bool CSG_Matrix::Add_Row(const double *Data)
1101 {
1102  if( Add_Rows(1) )
1103  {
1104  Set_Row(m_ny - 1, Data);
1105 
1106  return( true );
1107  }
1108 
1109  return( false );
1110 }
1111 
1113 {
1114  if( m_ny == 0 )
1115  {
1116  return( Create(Data.Get_Size(), 1) && Set_Row(0, Data.Get_Data()) );
1117  }
1118 
1119  return( m_nx <= Data.Get_Size() && Add_Row(Data.Get_Data()) );
1120 }
1121 
1122 //---------------------------------------------------------
1123 bool CSG_Matrix::Ins_Col(sLong Col, const double *Data)
1124 {
1125  if( Col >= 0 && Col <= m_nx )
1126  {
1127  CSG_Matrix Tmp(*this);
1128 
1129  if( Create(Tmp.m_nx + 1, Tmp.m_ny) )
1130  {
1131  for(sLong y=0; y<m_ny; y++)
1132  {
1133  double *pz = m_z[y], *pz_tmp = Tmp.m_z[y];
1134 
1135  for(sLong x=0; x<m_nx; x++, pz++)
1136  {
1137  if( x != Col )
1138  {
1139  *pz = *pz_tmp; pz_tmp++;
1140  }
1141  else if( Data )
1142  {
1143  *pz = Data[y];
1144  }
1145  }
1146  }
1147 
1148  return( true );
1149  }
1150  }
1151 
1152  return( false );
1153 }
1154 
1155 bool CSG_Matrix::Ins_Col(sLong Col, const CSG_Vector &Data)
1156 {
1157  return( m_nx == 0 ? Add_Col(Data) : m_ny <= Data.Get_Size() ? Ins_Col(Col, Data.Get_Data()) : false );
1158 }
1159 
1160 //---------------------------------------------------------
1161 bool CSG_Matrix::Ins_Row(sLong Row, const double *Data)
1162 {
1163  if( Row >= 0 && Row <= m_ny )
1164  {
1165  CSG_Matrix Tmp(*this);
1166 
1167  if( Create(Tmp.m_nx, Tmp.m_ny + 1) )
1168  {
1169  for(sLong y=0, y_tmp=0; y<m_ny; y++)
1170  {
1171  if( y != Row )
1172  {
1173  memcpy(m_z[y], Tmp.m_z[y_tmp++], m_nx * sizeof(double));
1174  }
1175  else if( Data )
1176  {
1177  memcpy(m_z[y], Data, m_nx * sizeof(double));
1178  }
1179  }
1180 
1181  return( true );
1182  }
1183  }
1184 
1185  return( false );
1186 }
1187 
1188 bool CSG_Matrix::Ins_Row(sLong Row, const CSG_Vector &Data)
1189 {
1190  return( m_ny == 0 ? Add_Row(Data) : m_nx <= Data.Get_Size() ? Ins_Row(Row, Data.Get_Data()) : false );
1191 }
1192 
1193 //---------------------------------------------------------
1194 bool CSG_Matrix::Set_Col(sLong Col, const double *Data)
1195 {
1196  if( Data && Col >= 0 && Col < m_nx )
1197  {
1198  for(int y=0; y<m_ny; y++)
1199  {
1200  m_z[y][Col] = Data[y];
1201  }
1202 
1203  return( true );
1204  }
1205 
1206  return( false );
1207 }
1208 
1209 bool CSG_Matrix::Set_Col(sLong Col, const CSG_Vector &Data)
1210 {
1211  return( m_ny <= Data.Get_Size() ? Set_Col(Col, Data.Get_Data()) : false );
1212 }
1213 
1214 //---------------------------------------------------------
1219 {
1220  return( Create(1, Data.Get_Size()) && Set_Col(0, Data.Get_Data()) );
1221 }
1222 
1223 //---------------------------------------------------------
1224 bool CSG_Matrix::Set_Row(sLong Row, const double *Data)
1225 {
1226  if( Data && Row >= 0 && Row < m_ny )
1227  {
1228  memcpy(m_z[Row], Data, m_nx * sizeof(double));
1229 
1230  return( true );
1231  }
1232 
1233  return( false );
1234 }
1235 
1236 bool CSG_Matrix::Set_Row(sLong Row, const CSG_Vector &Data)
1237 {
1238  return( m_nx <= Data.Get_Size() ? Set_Row(Row, Data.Get_Data()) : false );
1239 }
1240 
1241 //---------------------------------------------------------
1246 {
1247  return( Create(Data.Get_Size(), 1) && Set_Row(0, Data.Get_Data()) );
1248 }
1249 
1250 //---------------------------------------------------------
1252 {
1253  if( m_nx == 1 )
1254  {
1255  return( Destroy() );
1256  }
1257 
1258  if( Col >= 0 && Col < m_nx )
1259  {
1260  CSG_Matrix Tmp(*this);
1261 
1262  if( Create(Tmp.m_nx - 1, Tmp.m_ny) )
1263  {
1264  for(sLong y=0; y<m_ny; y++)
1265  {
1266  double *pz = m_z[y], *pz_tmp = Tmp.m_z[y];
1267 
1268  for(sLong x_tmp=0; x_tmp<Tmp.m_nx; x_tmp++, pz_tmp++)
1269  {
1270  if( x_tmp != Col )
1271  {
1272  *pz = *pz_tmp; pz++;
1273  }
1274  }
1275  }
1276 
1277  return( true );
1278  }
1279  }
1280 
1281  return( false );
1282 }
1283 
1284 //---------------------------------------------------------
1286 {
1287  if( m_ny == 1 )
1288  {
1289  return( Destroy() );
1290  }
1291 
1292  if( Row >= 0 && Row < m_ny )
1293  {
1294  CSG_Matrix Tmp(*this);
1295 
1296  if( Create(Tmp.m_nx, Tmp.m_ny - 1) )
1297  {
1298  for(sLong y=0, y_tmp=0; y_tmp<Tmp.m_ny; y_tmp++)
1299  {
1300  if( y_tmp != Row )
1301  {
1302  memcpy(m_z[y++], Tmp.m_z[y_tmp], m_nx * sizeof(double));
1303  }
1304  }
1305 
1306  return( true );
1307  }
1308  }
1309 
1310  return( false );
1311 }
1312 
1313 //---------------------------------------------------------
1315 {
1316  CSG_Vector Vector;
1317 
1318  if( Col >= 0 && Col < m_nx )
1319  {
1320  Vector.Create(m_ny);
1321 
1322  for(int y=0; y<m_ny; y++)
1323  {
1324  Vector[y] = m_z[y][Col];
1325  }
1326  }
1327 
1328  return( Vector );
1329 }
1330 
1331 //---------------------------------------------------------
1333 {
1334  CSG_Vector Vector;
1335 
1336  if( Row >= 0 && Row < m_ny )
1337  {
1338  Vector.Create(m_nx, m_z[Row]);
1339  }
1340 
1341  return( Vector );
1342 }
1343 
1344 
1346 // //
1348 
1349 //---------------------------------------------------------
1350 CSG_String CSG_Matrix::to_String(int Width, int Precision, bool bScientific, const SG_Char *Separator) const
1351 {
1352  CSG_String s, sep(Separator && *Separator ? Separator : SG_T(" "));
1353 
1354  int nDigits = SG_Get_Digit_Count(1 + (int)m_ny);
1355 
1356  for(sLong y=0; y<m_ny; y++)
1357  {
1358  s += CSG_String::Format("\n%0*ld:", nDigits, y + 1);
1359 
1360  for(sLong x=0; x<m_nx; x++)
1361  {
1362  s += sep + SG_Get_Double_asString(m_z[y][x], Width, Precision, bScientific);
1363  }
1364  }
1365 
1366  s += "\n";
1367 
1368  return( s );
1369 }
1370 
1371 //---------------------------------------------------------
1373 {
1374  Destroy();
1375 
1376  CSG_String_Tokenizer Lines(String, "\r\n");
1377 
1378  while( Lines.Has_More_Tokens() )
1379  {
1380  CSG_Vector Row; CSG_String_Tokenizer Line(Lines.Get_Next_Token().AfterFirst(':'));
1381 
1382  while( Line.Has_More_Tokens() )
1383  {
1384  double z; CSG_String s(Line.Get_Next_Token());
1385 
1386  if( s.asDouble(z) )
1387  {
1388  Row.Add_Row(z);
1389  }
1390  }
1391 
1392  Add_Row(Row);
1393  }
1394 
1395  return( Get_NRows() > 0 );
1396 }
1397 
1398 
1400 // //
1402 
1403 //---------------------------------------------------------
1404 bool CSG_Matrix::is_Equal(const CSG_Matrix &Matrix) const
1405 {
1406  if( m_nx < Matrix.m_nx || m_ny < Matrix.m_ny )
1407  {
1408  return( false );
1409  }
1410 
1411  for(sLong y=0; y<m_ny; y++) for(sLong x=0; x<m_nx; x++)
1412  {
1413  if( m_z[y][x] != Matrix.m_z[y][x] )
1414  {
1415  return( false );
1416  }
1417  }
1418 
1419  return( true );
1420 }
1421 
1422 //---------------------------------------------------------
1423 bool CSG_Matrix::Assign(double Scalar)
1424 {
1425  if( m_nx > 0 && m_ny > 0 )
1426  {
1427  for(sLong y=0; y<m_ny; y++)
1428  {
1429  for(sLong x=0; x<m_nx; x++)
1430  {
1431  m_z[y][x] = Scalar;
1432  }
1433  }
1434 
1435  return( true );
1436  }
1437 
1438  return( false );
1439 }
1440 
1441 //---------------------------------------------------------
1442 bool CSG_Matrix::Assign(const CSG_Matrix &Matrix)
1443 {
1444  return( Create(Matrix) );
1445 }
1446 
1447 //---------------------------------------------------------
1448 bool CSG_Matrix::Add(double Scalar)
1449 {
1450  if( m_nx > 0 && m_ny > 0 )
1451  {
1452  for(sLong y=0; y<m_ny; y++)
1453  {
1454  for(sLong x=0; x<m_nx; x++)
1455  {
1456  m_z[y][x] += Scalar;
1457  }
1458  }
1459 
1460  return( true );
1461  }
1462 
1463  return( false );
1464 }
1465 
1466 //---------------------------------------------------------
1467 bool CSG_Matrix::Add(const CSG_Matrix &Matrix)
1468 {
1469  if( m_nx == Matrix.m_nx && m_ny == Matrix.m_ny )
1470  {
1471  for(sLong y=0; y<m_ny; y++)
1472  {
1473  for(sLong x=0; x<m_nx; x++)
1474  {
1475  m_z[y][x] += Matrix.m_z[y][x];
1476  }
1477  }
1478 
1479  return( true );
1480  }
1481 
1482  return( false );
1483 }
1484 
1485 //---------------------------------------------------------
1487 {
1488  if( m_nx == Matrix.m_nx && m_ny == Matrix.m_ny )
1489  {
1490  for(sLong y=0; y<m_ny; y++)
1491  {
1492  for(sLong x=0; x<m_nx; x++)
1493  {
1494  m_z[y][x] -= Matrix.m_z[y][x];
1495  }
1496  }
1497 
1498  return( true );
1499  }
1500 
1501  return( false );
1502 }
1503 
1504 //---------------------------------------------------------
1505 bool CSG_Matrix::Multiply(double Scalar)
1506 {
1507  if( m_nx > 0 && m_ny > 0 )
1508  {
1509  for(sLong y=0; y<m_ny; y++)
1510  {
1511  for(sLong x=0; x<m_nx; x++)
1512  {
1513  m_z[y][x] *= Scalar;
1514  }
1515  }
1516 
1517  return( true );
1518  }
1519 
1520  return( false );
1521 }
1522 
1524 {
1525  CSG_Vector v;
1526 
1527  if( m_nx == Vector.Get_Size() && v.Create(m_ny) )
1528  {
1529  for(sLong y=0; y<m_ny; y++)
1530  {
1531  double z = 0.;
1532 
1533  for(sLong x=0; x<m_nx; x++)
1534  {
1535  z += m_z[y][x] * Vector(x);
1536  }
1537 
1538  v[y] = z;
1539  }
1540  }
1541 
1542  return( v );
1543 }
1544 
1546 {
1547  CSG_Matrix m;
1548 
1549  if( m_nx == Matrix.m_ny && m.Create(Matrix.m_nx, m_ny) )
1550  {
1551  for(int y=0; y<m.m_ny; y++)
1552  {
1553  for(int x=0; x<m.m_nx; x++)
1554  {
1555  double z = 0.;
1556 
1557  for(int n=0; n<m_nx; n++)
1558  {
1559  z += m_z[y][n] * Matrix.m_z[n][x];
1560  }
1561 
1562  m.m_z[y][x] = z;
1563  }
1564  }
1565  }
1566 
1567  return( m );
1568 }
1569 
1570 
1572 // //
1574 
1575 //---------------------------------------------------------
1576 bool CSG_Matrix::operator == (const CSG_Matrix &Matrix) const
1577 {
1578  return( is_Equal(Matrix) );
1579 }
1580 
1581 //---------------------------------------------------------
1583 {
1584  Assign(Scalar);
1585 
1586  return( *this );
1587 }
1588 
1590 {
1591  Assign(Matrix);
1592 
1593  return( *this );
1594 }
1595 
1596 //---------------------------------------------------------
1598 {
1599  Add(Scalar);
1600 
1601  return( *this );
1602 }
1603 
1605 {
1606  Add(Matrix);
1607 
1608  return( *this );
1609 }
1610 
1611 //---------------------------------------------------------
1613 {
1614  Add(-Scalar);
1615 
1616  return( *this );
1617 }
1618 
1620 {
1621  Subtract(Matrix);
1622 
1623  return( *this );
1624 }
1625 
1626 //---------------------------------------------------------
1628 {
1629  Multiply(Scalar);
1630 
1631  return( *this );
1632 }
1633 
1635 {
1636  Multiply(Matrix);
1637 
1638  return( *this );
1639 }
1640 
1641 //---------------------------------------------------------
1643 {
1644  CSG_Matrix m(*this);
1645 
1646  m.Add(Scalar);
1647 
1648  return( m );
1649 }
1650 
1652 {
1653  CSG_Matrix m(*this);
1654 
1655  m.Add(Matrix);
1656 
1657  return( m );
1658 }
1659 
1660 //---------------------------------------------------------
1662 {
1663  CSG_Matrix m(*this);
1664 
1665  m.Add(-Scalar);
1666 
1667  return( m );
1668 }
1669 
1671 {
1672  CSG_Matrix m(*this);
1673 
1674  m.Subtract(Matrix);
1675 
1676  return( m );
1677 }
1678 
1679 //---------------------------------------------------------
1681 {
1682  CSG_Matrix m(*this);
1683 
1684  m.Multiply(Scalar);
1685 
1686  return( m );
1687 }
1688 
1690 {
1691  return( Multiply(Vector) );
1692 }
1693 
1695 {
1696  return( Multiply(Matrix) );
1697 }
1698 
1699 CSG_Matrix operator * (double Scalar, const CSG_Matrix &Matrix)
1700 {
1701  return( Matrix * Scalar );
1702 }
1703 
1704 
1706 // //
1708 
1709 //---------------------------------------------------------
1711 {
1712  return( Create(m_nx, m_ny) );
1713 }
1714 
1715 //---------------------------------------------------------
1717 {
1718  if( m_nx > 0 && m_ny > 0 )
1719  {
1720  for(int y=0; y<m_ny; y++)
1721  {
1722  for(int x=0; x<m_nx; x++)
1723  {
1724  m_z[y][x] = x == y ? 1. : 0.;
1725  }
1726  }
1727 
1728  return( true );
1729  }
1730 
1731  return( false );
1732 }
1733 
1734 //---------------------------------------------------------
1736 {
1737  CSG_Matrix m;
1738 
1739  if( m.Create(*this) && Create(m_ny, m_nx) )
1740  {
1741  for(int y=0; y<m_ny; y++)
1742  {
1743  for(int x=0; x<m_nx; x++)
1744  {
1745  m_z[y][x] = m.m_z[x][y];
1746  }
1747  }
1748 
1749  return( true );
1750  }
1751 
1752  return( false );
1753 }
1754 
1755 //---------------------------------------------------------
1756 bool CSG_Matrix::Set_Inverse(bool bSilent, int nSubSquare)
1757 {
1758  int n = 0;
1759 
1760  if( nSubSquare > 0 )
1761  {
1762  if( nSubSquare <= m_nx && nSubSquare <= m_ny )
1763  {
1764  n = nSubSquare;
1765  }
1766  }
1767  else if( is_Square() )
1768  {
1769  n = (int)m_nx;
1770  }
1771 
1772  //-----------------------------------------------------
1773  if( n > 0 )
1774  {
1775  CSG_Matrix m(*this); CSG_Array p(sizeof(int), n);
1776 
1777  if( SG_Matrix_LU_Decomposition(n, (int *)p.Get_Array(), m.Get_Data(), bSilent) )
1778  {
1779  CSG_Vector v(n);
1780 
1781  for(int j=0; j<n && (bSilent || SG_UI_Process_Set_Progress(j, n)); j++)
1782  {
1783  v.Set_Zero(); v[j] = 1.;
1784 
1785  SG_Matrix_LU_Solve(n, (int *)p.Get_Array(), m, v.Get_Data(), true);
1786 
1787  for(int i=0; i<n; i++)
1788  {
1789  m_z[i][j] = v[i];
1790  }
1791  }
1792 
1793  return( true );
1794  }
1795  }
1796 
1797  return( false );
1798 }
1799 
1800 
1802 // //
1804 
1805 //---------------------------------------------------------
1807 {
1808  double d = 0.;
1809 
1810  if( is_Square() ) // det is only defined for squared matrices !
1811  {
1812  int n; CSG_Matrix m(*this); CSG_Array p(sizeof(int), Get_NX());
1813 
1814  if( SG_Matrix_LU_Decomposition(Get_NX(), (int *)p.Get_Array(), m.Get_Data(), true, &n) )
1815  {
1816  d = n % 2 ? -1. : 1.;
1817 
1818  for(int i=0; i<Get_NX(); i++)
1819  {
1820  d *= m[i][i];
1821  }
1822  }
1823  }
1824 
1825  return( d );
1826 }
1827 
1828 //---------------------------------------------------------
1830 {
1831  CSG_Matrix m(m_ny, m_nx);
1832 
1833  for(sLong y=0; y<m_ny; y++)
1834  {
1835  for(sLong x=0; x<m_nx; x++)
1836  {
1837  m.m_z[x][y] = m_z[y][x];
1838  }
1839  }
1840 
1841  return( m );
1842 }
1843 
1844 //---------------------------------------------------------
1845 CSG_Matrix CSG_Matrix::Get_Inverse(bool bSilent, int nSubSquare) const
1846 {
1847  CSG_Matrix m(*this);
1848 
1849  m.Set_Inverse(bSilent, nSubSquare);
1850 
1851  return( m );
1852 }
1853 
1854 
1856 // //
1857 // //
1858 // //
1860 
1861 //---------------------------------------------------------
1866 CSG_Matrix SG_Matrix_Get_Rotation(double R, bool bDegree)
1867 {
1868  if( bDegree ) { R *= M_DEG_TO_RAD; }
1869 
1870  double sin_a = sin(R), cos_a = cos(R);
1871 
1872  CSG_Matrix M(2, 2);
1873 
1874  M[0][0] = cos_a; M[0][1] = -sin_a;
1875  M[1][0] = sin_a; M[1][1] = cos_a;
1876 
1877  return( M );
1878 }
1879 
1880 //---------------------------------------------------------
1885 CSG_Matrix SG_Matrix_Get_Rotation(double Rx, double Ry, double Rz, bool bDegree)
1886 {
1887  if( bDegree ) { Rx *= M_DEG_TO_RAD; Ry *= M_DEG_TO_RAD; Rz *= M_DEG_TO_RAD; }
1888 
1889  double sin_a = sin(Rx), cos_a = cos(Rx);
1890  double sin_b = sin(Ry), cos_b = cos(Ry);
1891  double sin_c = sin(Rz), cos_c = cos(Rz);
1892 
1893  CSG_Matrix M(3, 3);
1894 
1895  M[0][0] = cos_a * cos_b; M[0][1] = cos_a * sin_b * sin_c - sin_a * cos_c; M[0][2] = cos_a * sin_b * cos_c + sin_a * sin_c;
1896  M[1][0] = sin_a * cos_b; M[1][1] = sin_a * sin_b * sin_c + cos_a * cos_c; M[1][2] = sin_a * sin_b * cos_c - cos_a * sin_c;
1897  M[2][0] = -sin_b; M[2][1] = cos_b * sin_c ; M[2][2] = cos_b * cos_c ;
1898 
1899  return( M );
1900 }
1901 
1902 
1904 // //
1905 // //
1906 // //
1908 
1909 //---------------------------------------------------------
1910 bool SG_Matrix_Solve(CSG_Matrix &Matrix, CSG_Vector &Vector, bool bSilent)
1911 {
1912  int n = Vector.Get_N(); // ...it makes no sense to perform this operation on long long arrays! ...maybe in (far) future ;-)
1913 
1914  if( n > 0 && n == Matrix.Get_NX() && n == Matrix.Get_NY() )
1915  {
1916  CSG_Array Permutation(sizeof(int), n);
1917 
1918  if( SG_Matrix_LU_Decomposition(n, (int *)Permutation.Get_Array(), Matrix.Get_Data(), bSilent) )
1919  {
1920  return( SG_Matrix_LU_Solve(n, (int *)Permutation.Get_Array(), Matrix, Vector.Get_Data(), bSilent) );
1921  }
1922  }
1923 
1924  return( false );
1925 }
1926 
1927 
1929 // //
1931 
1932 //---------------------------------------------------------
1933 bool SG_Matrix_Eigen_Reduction(const CSG_Matrix &Matrix, CSG_Matrix &Eigen_Vectors, CSG_Vector &Eigen_Values, bool bSilent)
1934 {
1935  CSG_Vector Intermediate; Eigen_Vectors = Matrix;
1936 
1937  return( SG_Matrix_Triangular_Decomposition(Eigen_Vectors, Eigen_Values, Intermediate) // Triangular decomposition (Householder's method)
1938  && SG_Matrix_Tridiagonal_QL (Eigen_Vectors, Eigen_Values, Intermediate) // Reduction of symmetric tridiagonal matrix
1939  );
1940 }
1941 
1942 
1944 // //
1946 
1947 //---------------------------------------------------------
1948 bool SG_Matrix_LU_Decomposition(int n, int *Permutation, double **Matrix, bool bSilent, int *nRowChanges)
1949 {
1950  CSG_Vector Vector(n); if( nRowChanges ) { (*nRowChanges) = 0; }
1951 
1952  for(int i=0; i<n && (bSilent || SG_UI_Process_Set_Progress(i, n)); i++)
1953  {
1954  double dMax = 0.;
1955 
1956  for(int j=0; j<n; j++)
1957  {
1958  double d = fabs(Matrix[i][j]);
1959 
1960  if( d > dMax )
1961  {
1962  dMax = d;
1963  }
1964  }
1965 
1966  if( dMax <= 0. ) // singular matrix !!!...
1967  {
1968  return( false );
1969  }
1970 
1971  Vector[i] = 1. / dMax;
1972  }
1973 
1974  int iMax = 0;
1975 
1976  for(int j=0; j<n && (bSilent || SG_UI_Process_Set_Progress(j, n)); j++)
1977  {
1978  for(int i=0; i<j; i++)
1979  {
1980  double Sum = Matrix[i][j];
1981 
1982  for(int k=0; k<i; k++)
1983  {
1984  Sum -= Matrix[i][k] * Matrix[k][j];
1985  }
1986 
1987  Matrix[i][j] = Sum;
1988  }
1989 
1990  double dMax = 0.;
1991 
1992  for(int i=j; i<n; i++)
1993  {
1994  double Sum = Matrix[i][j];
1995 
1996  for(int k=0; k<j; k++)
1997  {
1998  Sum -= Matrix[i][k] * Matrix[k][j];
1999  }
2000 
2001  Matrix[i][j] = Sum;
2002 
2003  double d = Vector[i] * fabs(Sum);
2004 
2005  if( d >= dMax )
2006  {
2007  dMax = d;
2008  iMax = i;
2009  }
2010  }
2011 
2012  if( j != iMax )
2013  {
2014  for(int k=0; k<n; k++)
2015  {
2016  double d = Matrix[iMax][k];
2017  Matrix[iMax][k] = Matrix[j ][k];
2018  Matrix[j ][k] = d;
2019  }
2020 
2021  Vector[iMax] = Vector[j];
2022 
2023  if( nRowChanges ) { (*nRowChanges)++; }
2024  }
2025 
2026  Permutation[j] = iMax;
2027 
2028  if( Matrix[j][j] == 0. )
2029  {
2030  Matrix[j][j] = M_FLT_EPSILON;
2031  }
2032 
2033  if( j != n )
2034  {
2035  double d = 1. / (Matrix[j][j]);
2036 
2037  for(int i=j+1; i<n; i++)
2038  {
2039  Matrix[i][j] *= d;
2040  }
2041  }
2042  }
2043 
2044  return( bSilent || SG_UI_Process_Get_Okay(false) );
2045 }
2046 
2047 
2049 // //
2051 
2052 //---------------------------------------------------------
2053 bool SG_Matrix_LU_Solve(int n, const int *Permutation, const double **Matrix, double *Vector, bool bSilent)
2054 {
2055  for(int i=0, k=-1; i<n && (bSilent || SG_UI_Process_Set_Progress(i, n)); i++)
2056  {
2057  double Sum = Vector[Permutation[i]]; Vector[Permutation[i]] = Vector[i];
2058 
2059  if( k >= 0 )
2060  {
2061  for(int j=k; j<=i-1; j++)
2062  {
2063  Sum -= Matrix[i][j] * Vector[j];
2064  }
2065  }
2066  else if( Sum )
2067  {
2068  k = i;
2069  }
2070 
2071  Vector[i] = Sum;
2072  }
2073 
2074  for(int i=n-1; i>=0 && (bSilent || SG_UI_Process_Set_Progress(n-i, n)); i--)
2075  {
2076  double Sum = Vector[i];
2077 
2078  for(int j=i+1; j<n; j++)
2079  {
2080  Sum -= Matrix[i][j] * Vector[j];
2081  }
2082 
2083  Vector[i] = Sum / Matrix[i][i];
2084  }
2085 
2086  return( true );
2087 }
2088 
2089 
2091 // //
2093 
2094 //---------------------------------------------------------
2095 // Householder reduction of matrix a to tridiagonal form.
2096 
2098 {
2099  if( A.Get_NX() != A.Get_NY() )
2100  {
2101  return( false );
2102  }
2103 
2104  int n = A.Get_NX();
2105 
2106  if( !d.Create(n) || !e.Create(n) )
2107  {
2108  return( false );
2109  }
2110 
2111  for(int i=n-1, l=n-2; i>=1; i--, l--)
2112  {
2113  double h = 0.;
2114 
2115  if( l > 0 )
2116  {
2117  double scale = 0.;
2118 
2119  for(int k=0; k<=l; k++)
2120  {
2121  scale += fabs(A[i][k]);
2122  }
2123 
2124  if( scale == 0. )
2125  {
2126  e[i] = A[i][l];
2127  }
2128  else
2129  {
2130  for(int k=0; k<=l; k++)
2131  {
2132  A[i][k] /= scale;
2133  h += A[i][k] * A[i][k];
2134  }
2135 
2136  double f = A[i][l];
2137  double g = f > 0. ? -sqrt(h) : sqrt(h);
2138  e[i] = scale * g;
2139  h -= f * g;
2140  A[i][l] = f - g;
2141  f = 0.;
2142 
2143  for(int j=0; j<=l; j++)
2144  {
2145  A[j][i] = A[i][j]/h;
2146  g = 0.;
2147 
2148  for(int k=0; k<=j; k++)
2149  {
2150  g += A[j][k] * A[i][k];
2151  }
2152 
2153  for(int k=j+1; k<=l; k++)
2154  {
2155  g += A[k][j] * A[i][k];
2156  }
2157 
2158  e[j] = g / h;
2159  f += e[j] * A[i][j];
2160  }
2161 
2162  double hh = f / (h + h);
2163 
2164  for(int j=0; j<=l; j++)
2165  {
2166  f = A[i][j];
2167  e[j] = g = e[j] - hh * f;
2168 
2169  for(int k=0; k<=j; k++)
2170  {
2171  A[j][k] -= (f * e[k] + g * A[i][k]);
2172  }
2173  }
2174  }
2175  }
2176  else
2177  {
2178  e[i] = A[i][l];
2179  }
2180 
2181  d[i] = h;
2182  }
2183 
2184  d[0] = 0.;
2185  e[0] = 0.;
2186 
2187  for(int i=0, l=-1; i<n; i++, l++)
2188  {
2189  if( d[i] )
2190  {
2191  for(int j=0; j<=l; j++)
2192  {
2193  double g = 0.;
2194 
2195  for(int k=0; k<=l; k++)
2196  {
2197  g += A[i][k] * A[k][j];
2198  }
2199 
2200  for(int k=0; k<=l; k++)
2201  {
2202  A[k][j] -= g * A[k][i];
2203  }
2204  }
2205  }
2206 
2207  d[i] = A[i][i];
2208  A[i][i] = 1.;
2209 
2210  for(int j=0; j<=l; j++)
2211  {
2212  A[j][i] = A[i][j] = 0.;
2213  }
2214  }
2215 
2216  return( true );
2217 }
2218 
2219 
2221 // //
2223 
2224 //---------------------------------------------------------
2225 // Tridiagonal QL algorithm -- Implicit
2226 
2228 {
2229  if( Q.Get_NX() != Q.Get_NY() || Q.Get_NX() != d.Get_N() || Q.Get_NX() != e.Get_N() )
2230  {
2231  return( false );
2232  }
2233 
2234  int m, l, iter, i, k, n;
2235  double s, r, p, g, f, dd, c, b;
2236 
2237  n = d.Get_N();
2238 
2239  for(i=1; i<n; i++)
2240  {
2241  e[i - 1] = e[i];
2242  }
2243 
2244  e[n - 1] = 0.;
2245 
2246  for(l=0; l<n; l++)
2247  {
2248  iter = 0;
2249 
2250  do
2251  {
2252  for(m=l; m<n-1; m++)
2253  {
2254  dd = fabs(d[m]) + fabs(d[m + 1]);
2255 
2256  if( fabs(e[m]) + dd == dd )
2257  {
2258  break;
2259  }
2260  }
2261 
2262  if( m != l )
2263  {
2264  if( iter++ == 30 )
2265  {
2266  return( false ); // no convergence in TLQI !
2267  }
2268 
2269  g = (d[l+1] - d[l]) / (2. * e[l]);
2270  r = sqrt((g * g) + 1.);
2271  g = d[m] - d[l] + e[l] / (g + M_SET_SIGN(r, g));
2272  s = c = 1.;
2273  p = 0.;
2274 
2275  for(i = m-1; i >= l; i--)
2276  {
2277  f = s * e[i];
2278  b = c * e[i];
2279 
2280  if (fabs(f) >= fabs(g))
2281  {
2282  c = g / f;
2283  r = sqrt((c * c) + 1.);
2284  e[i+1] = f * r;
2285  c *= (s = 1. / r);
2286  }
2287  else
2288  {
2289  s = f / g;
2290  r = sqrt((s * s) + 1.);
2291  e[i+1] = g * r;
2292  s *= (c = 1. / r);
2293  }
2294 
2295  g = d[i+1] - p;
2296  r = (d[i] - g) * s + 2. * c * b;
2297  p = s * r;
2298  d[i+1] = g + p;
2299  g = c * r - b;
2300 
2301  for(k=0; k<n; k++)
2302  {
2303  f = Q[k][i+1];
2304  Q[k][i+1] = s * Q[k][i] + c * f;
2305  Q[k][i] = c * Q[k][i] - s * f;
2306  }
2307  }
2308 
2309  d[l] = d[l] - p;
2310  e[l] = g;
2311  e[m] = 0.;
2312  }
2313  }
2314  while( m != l );
2315  }
2316 
2317  return( true );
2318 }
2319 
2320 
2322 // //
2323 // //
2324 // //
2326 
2327 //---------------------------------------------------------
M_DEG_TO_RAD
#define M_DEG_TO_RAD
Definition: mat_tools.h:109
CSG_Array::Set_Array
bool Set_Array(sLong nValues, bool bShrink=true)
Definition: api_memory.cpp:310
CSG_Matrix::to_String
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
Definition: mat_matrix.cpp:1350
CSG_Vector::from_String
bool from_String(const CSG_String &String)
Definition: mat_matrix.cpp:238
CSG_Matrix::Get_Inverse
CSG_Matrix Get_Inverse(bool bSilent=true, int nSubSquare=0) const
Definition: mat_matrix.cpp:1845
CSG_Matrix::is_Equal
bool is_Equal(const CSG_Matrix &Matrix) const
Definition: mat_matrix.cpp:1404
SG_Get_Double_asString
CSG_String SG_Get_Double_asString(double Number, int Width, int Precision, bool bScientific)
Definition: mat_tools.cpp:159
CSG_Vector::is_Collinear
bool is_Collinear(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:296
SG_T
#define SG_T(s)
Definition: api_core.h:537
CSG_String_Tokenizer::Has_More_Tokens
bool Has_More_Tokens(void) const
Definition: api_string.cpp:1511
CSG_Matrix::Set_Identity
bool Set_Identity(void)
Definition: mat_matrix.cpp:1716
CSG_Vector::Get_Data
double * Get_Data(void) const
Definition: mat_tools.h:384
CSG_Vector::operator=
CSG_Vector & operator=(double Scalar)
Definition: mat_matrix.cpp:460
CSG_Matrix::Del_Rows
bool Del_Rows(sLong nRows)
Definition: mat_matrix.cpp:1061
CSG_Vector::to_String
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
Definition: mat_matrix.cpp:225
CSG_Matrix::Set_Inverse
bool Set_Inverse(bool bSilent=true, int nSubSquare=0)
Definition: mat_matrix.cpp:1756
M_SET_SIGN
#define M_SET_SIGN(x, sign)
Definition: mat_tools.h:145
SG_VectorR2_Rotate
bool SG_VectorR2_Rotate(double &x, double &y, double Angle)
Definition: mat_matrix.cpp:750
A
#define A
CSG_Matrix::Get_NRows
sLong Get_NRows(void) const
Definition: mat_tools.h:523
CSG_Matrix::Del_Col
bool Del_Col(sLong Col)
Definition: mat_matrix.cpp:1251
CSG_Matrix::~CSG_Matrix
virtual ~CSG_Matrix(void)
Definition: mat_matrix.cpp:918
CSG_Matrix::Add
bool Add(double Scalar)
Definition: mat_matrix.cpp:1448
CSG_Matrix::CSG_Matrix
CSG_Matrix(void)
Definition: mat_matrix.cpp:823
CSG_Vector::Get_Unity
CSG_Vector Get_Unity(void) const
Definition: mat_matrix.cpp:733
CSG_Vector::Del_Row
bool Del_Row(sLong Row=-1)
Definition: mat_matrix.cpp:201
SG_Malloc
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
Definition: api_memory.cpp:65
CSG_Vector::Set_Unity
bool Set_Unity(void)
Definition: mat_matrix.cpp:596
SG_UI_Process_Get_Okay
bool SG_UI_Process_Get_Okay(bool bBlink)
Definition: api_callback.cpp:207
CSG_Matrix::operator+
CSG_Matrix operator+(double Scalar) const
Definition: mat_matrix.cpp:1642
CSG_Matrix::Set_Cols
bool Set_Cols(sLong nCols)
Definition: mat_matrix.cpp:958
CSG_Matrix::Add_Rows
bool Add_Rows(sLong nRows)
Definition: mat_matrix.cpp:1011
CSG_Vector::Add_Row
bool Add_Row(double Value=0.)
Definition: mat_matrix.cpp:188
CSG_Vector::Get_Cross_Product
CSG_Vector Get_Cross_Product(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:675
SG_VectorR3_Rotate
bool SG_VectorR3_Rotate(double Vector[3], size_t Axis, double Angle)
Definition: mat_matrix.cpp:776
SG_Free
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
Definition: api_memory.cpp:83
CSG_Vector::operator*
CSG_Vector operator*(double Scalar) const
Definition: mat_matrix.cpp:565
CSG_Matrix::Del_Row
bool Del_Row(sLong Row)
Definition: mat_matrix.cpp:1285
CSG_Vector::operator-
CSG_Vector operator-(double Scalar) const
Definition: mat_matrix.cpp:546
SG_Compare_Double
int SG_Compare_Double(const void *a, const void *b)
Definition: mat_tools.cpp:199
CSG_Matrix::Get_Col
CSG_Vector Get_Col(sLong Col) const
Definition: mat_matrix.cpp:1314
CSG_Matrix::Set_Row
bool Set_Row(sLong Row, const double *Data)
Definition: mat_matrix.cpp:1224
CSG_Vector::Get_Length
double Get_Length(void) const
Definition: mat_matrix.cpp:685
CSG_Vector::CSG_Vector
CSG_Vector(void)
Definition: mat_matrix.cpp:74
CSG_Matrix::from_String
bool from_String(const CSG_String &String)
Definition: mat_matrix.cpp:1372
CSG_Matrix::Get_Data
double ** Get_Data(void) const
Definition: mat_tools.h:526
CSG_Matrix::Add_Cols
bool Add_Cols(sLong nCols)
Definition: mat_matrix.cpp:990
mat_tools.h
CSG_Matrix::Create
bool Create(const CSG_Matrix &Matrix)
Definition: mat_matrix.cpp:836
CSG_Vector::Flip_Values
bool Flip_Values(void)
Definition: mat_matrix.cpp:619
CSG_Vector
Definition: mat_tools.h:360
CSG_Vector::Set_Zero
bool Set_Zero(void)
Definition: mat_matrix.cpp:590
CSG_Matrix::Set_Transpose
bool Set_Transpose(void)
Definition: mat_matrix.cpp:1735
CSG_Matrix::Multiply
bool Multiply(double Scalar)
Definition: mat_matrix.cpp:1505
CSG_Vector::operator+
CSG_Vector operator+(double Scalar) const
Definition: mat_matrix.cpp:527
SG_Matrix_Triangular_Decomposition
bool SG_Matrix_Triangular_Decomposition(CSG_Matrix &A, CSG_Vector &d, CSG_Vector &e)
Definition: mat_matrix.cpp:2097
CSG_Vector::Create
bool Create(const CSG_Vector &Vector)
Definition: mat_matrix.cpp:87
CSG_Vector::Assign
bool Assign(double Scalar)
Definition: mat_matrix.cpp:317
sLong
signed long long sLong
Definition: api_core.h:158
CSG_Matrix::Assign
bool Assign(double Scalar)
Definition: mat_matrix.cpp:1423
CSG_Matrix::Subtract
bool Subtract(const CSG_Matrix &Matrix)
Definition: mat_matrix.cpp:1486
SG_Matrix_LU_Decomposition
bool SG_Matrix_LU_Decomposition(int n, int *Permutation, double **Matrix, bool bSilent, int *nRowChanges)
Definition: mat_matrix.cpp:1948
CSG_Vector::Get_Angle
double Get_Angle(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:703
CSG_Vector::Get_Scalar_Product
CSG_Vector Get_Scalar_Product(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:665
CSG_Vector::Get_N
int Get_N(void) const
Definition: mat_tools.h:382
CSG_Vector::operator-=
CSG_Vector & operator-=(double Scalar)
Definition: mat_matrix.cpp:490
CSG_Vector::Multiply_Cross
bool Multiply_Cross(const CSG_Vector &Vector)
Definition: mat_matrix.cpp:416
CSG_Matrix::Del_Cols
bool Del_Cols(sLong nCols)
Definition: mat_matrix.cpp:1037
CSG_Vector::Subtract
bool Subtract(const CSG_Vector &Vector)
Definition: mat_matrix.cpp:378
SG_Matrix_LU_Solve
bool SG_Matrix_LU_Solve(int n, const int *Permutation, const double **Matrix, double *Vector, bool bSilent)
Definition: mat_matrix.cpp:2053
CSG_Array::Create
void * Create(const CSG_Array &Array)
Definition: api_memory.cpp:250
CSG_Matrix::Set_Rows
bool Set_Rows(sLong nRows)
Definition: mat_matrix.cpp:974
CSG_Vector::is_Null
bool is_Null(void) const
Definition: mat_matrix.cpp:263
M_FLT_EPSILON
#define M_FLT_EPSILON
Definition: mat_tools.h:121
CSG_Vector::~CSG_Vector
virtual ~CSG_Vector(void)
Definition: mat_matrix.cpp:125
CSG_Array::Get_Array
void * Get_Array(void) const
Definition: api_core.h:336
CSG_Matrix::Set_Col
bool Set_Col(sLong Col, const double *Data)
Definition: mat_matrix.cpp:1194
SG_Matrix_Solve
bool SG_Matrix_Solve(CSG_Matrix &Matrix, CSG_Vector &Vector, bool bSilent)
Definition: mat_matrix.cpp:1910
CSG_String::Format
static CSG_String Format(const char *Format,...)
Definition: api_string.cpp:270
SG_Matrix_Eigen_Reduction
bool SG_Matrix_Eigen_Reduction(const CSG_Matrix &Matrix, CSG_Matrix &Eigen_Vectors, CSG_Vector &Eigen_Values, bool bSilent)
Definition: mat_matrix.cpp:1933
CSG_Matrix::Get_Transpose
CSG_Matrix Get_Transpose(void) const
Definition: mat_matrix.cpp:1829
CSG_Vector::Multiply
bool Multiply(double Scalar)
Definition: mat_matrix.cpp:394
CSG_Matrix::is_Square
bool is_Square(void) const
Definition: mat_tools.h:536
CSG_Vector::Get_Size
sLong Get_Size(void) const
Definition: mat_tools.h:380
CSG_String::AfterFirst
CSG_String AfterFirst(char Character) const
Definition: api_string.cpp:644
CSG_Vector::is_Equal
bool is_Equal(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:277
CSG_Vector::Destroy
bool Destroy(void)
Definition: mat_matrix.cpp:130
CSG_Vector::Add
bool Add(double Scalar)
Definition: mat_matrix.cpp:346
SG_Char
#define SG_Char
Definition: api_core.h:536
CSG_Array::Inc_Array
bool Inc_Array(sLong nValues=1)
Definition: api_memory.cpp:414
CSG_Array
Definition: api_core.h:308
CSG_Matrix::operator+=
CSG_Matrix & operator+=(double Scalar)
Definition: mat_matrix.cpp:1597
B
#define B
CSG_String
Definition: api_core.h:563
SG_Get_Digit_Count
int SG_Get_Digit_Count(int Number)
Definition: mat_tools.cpp:144
CSG_Array::Dec_Array
bool Dec_Array(bool bShrink=true)
Definition: api_memory.cpp:425
SG_Matrix_Get_Rotation
CSG_Matrix SG_Matrix_Get_Rotation(double R, bool bDegree)
Definition: mat_matrix.cpp:1866
CSG_Matrix::Set_Zero
bool Set_Zero(void)
Definition: mat_matrix.cpp:1710
CSG_Matrix::Ins_Col
bool Ins_Col(sLong Col, const double *Data=NULL)
Definition: mat_matrix.cpp:1123
SG_UI_Process_Set_Progress
bool SG_UI_Process_Set_Progress(int Position, int Range)
Definition: api_callback.cpp:255
CSG_Vector::Del_Rows
bool Del_Rows(sLong nRows)
Definition: mat_matrix.cpp:177
CSG_Matrix::operator-
CSG_Matrix operator-(double Scalar) const
Definition: mat_matrix.cpp:1661
CSG_Matrix::Destroy
bool Destroy(void)
Definition: mat_matrix.cpp:923
CSG_Matrix::Get_Row
CSG_Vector Get_Row(sLong Row) const
Definition: mat_matrix.cpp:1332
CSG_Matrix::operator==
bool operator==(const CSG_Matrix &Matrix) const
Definition: mat_matrix.cpp:1576
CSG_Matrix::Add_Col
bool Add_Col(const double *Data=NULL)
Definition: mat_matrix.cpp:1077
CSG_Matrix::operator*=
CSG_Matrix & operator*=(double Scalar)
Definition: mat_matrix.cpp:1627
CSG_Matrix::Get_Determinant
double Get_Determinant(void) const
Definition: mat_matrix.cpp:1806
CSG_Matrix::Add_Row
bool Add_Row(const double *Data=NULL)
Definition: mat_matrix.cpp:1100
CSG_Vector::operator+=
CSG_Vector & operator+=(double Scalar)
Definition: mat_matrix.cpp:475
CSG_Matrix::operator=
CSG_Matrix & operator=(double Scalar)
Definition: mat_matrix.cpp:1582
SG_Matrix_Tridiagonal_QL
bool SG_Matrix_Tridiagonal_QL(CSG_Matrix &Q, CSG_Vector &d, CSG_Vector &e)
Definition: mat_matrix.cpp:2227
CSG_String::asDouble
double asDouble(void) const
Definition: api_string.cpp:760
CSG_Vector::Multiply_Scalar
double Multiply_Scalar(const CSG_Vector &Vector) const
Definition: mat_matrix.cpp:433
SG_Realloc
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
Definition: api_memory.cpp:77
CSG_Vector::operator*=
CSG_Vector & operator*=(double Scalar)
Definition: mat_matrix.cpp:505
CSG_Matrix::Set_Size
bool Set_Size(sLong nRows, sLong nCols)
Definition: mat_matrix.cpp:952
operator*
CSG_Vector operator*(double Scalar, const CSG_Vector &Vector)
Definition: mat_matrix.cpp:579
CSG_Vector::Set_Rows
bool Set_Rows(sLong nRows)
Definition: mat_matrix.cpp:140
CSG_String_Tokenizer
Definition: api_core.h:760
CSG_Matrix
Definition: mat_tools.h:478
CSG_Matrix::operator*
CSG_Matrix operator*(double Scalar) const
Definition: mat_matrix.cpp:1680
CSG_String_Tokenizer::Get_Next_Token
CSG_String Get_Next_Token(void)
Definition: api_string.cpp:1489
CSG_Vector::Sort
bool Sort(bool bAscending=true)
Definition: mat_matrix.cpp:642
CSG_Vector::Add_Rows
bool Add_Rows(sLong nRows)
Definition: mat_matrix.cpp:156
CSG_Matrix::Get_NY
int Get_NY(void) const
Definition: mat_tools.h:521
CSG_Matrix::Ins_Row
bool Ins_Row(sLong Row, const double *Data=NULL)
Definition: mat_matrix.cpp:1161
CSG_Matrix::Get_NX
int Get_NX(void) const
Definition: mat_tools.h:520
CSG_Matrix::operator-=
CSG_Matrix & operator-=(double Scalar)
Definition: mat_matrix.cpp:1612