SAGA API v9.10
Loading...
Searching...
No Matches
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
87bool CSG_Vector::Create(const CSG_Vector &Vector)
88{
89 return( Assign(Vector) );
90}
91
92//---------------------------------------------------------
93CSG_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
100bool 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//---------------------------------------------------------
188bool 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//---------------------------------------------------------
225CSG_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//---------------------------------------------------------
263bool 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//---------------------------------------------------------
277bool 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//---------------------------------------------------------
296bool 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//---------------------------------------------------------
317bool 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//---------------------------------------------------------
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//---------------------------------------------------------
346bool 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//---------------------------------------------------------
362bool 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//---------------------------------------------------------
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//---------------------------------------------------------
394bool 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//---------------------------------------------------------
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//---------------------------------------------------------
433double 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//---------------------------------------------------------
449bool 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
574double CSG_Vector::operator * (const CSG_Vector &Vector) const
575{
576 return( Multiply_Scalar(Vector) );
577}
578
579CSG_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//---------------------------------------------------------
642bool 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//---------------------------------------------------------
685double 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//---------------------------------------------------------
703double 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//---------------------------------------------------------
750bool 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//---------------------------------------------------------
764bool SG_VectorR2_Rotate(double Vector[2], double Angle)
765{
766 return( SG_VectorR2_Rotate(Vector[0], Vector[1], Angle) );
767}
768
769//---------------------------------------------------------
770bool 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//---------------------------------------------------------
776bool 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//---------------------------------------------------------
810bool 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
837{
838 return( Create(Matrix.m_nx, Matrix.m_ny, Matrix.m_z[0]) );
839}
840
841//---------------------------------------------------------
842CSG_Matrix::CSG_Matrix(sLong nCols, sLong nRows, const double *Data)
843{
844 _On_Construction();
845
846 Create(nCols, nRows, Data);
847}
848
849bool 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//---------------------------------------------------------
892CSG_Matrix::CSG_Matrix(sLong nCols, sLong nRows, const double **Data)
893{
894 _On_Construction();
895
896 Create(nCols, nRows, Data);
897}
898
899bool 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//---------------------------------------------------------
941void CSG_Matrix::_On_Construction(void)
942{
943 m_z = NULL; m_nx = m_ny = 0;
944}
945
946
948// //
950
951//---------------------------------------------------------
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//---------------------------------------------------------
1077bool 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//---------------------------------------------------------
1100bool 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//---------------------------------------------------------
1123bool 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
1156{
1157 return( m_nx == 0 ? Add_Col(Data) : m_ny <= Data.Get_Size() ? Ins_Col(Col, Data.Get_Data()) : false );
1158}
1159
1160//---------------------------------------------------------
1161bool 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
1189{
1190 return( m_ny == 0 ? Add_Row(Data) : m_nx <= Data.Get_Size() ? Ins_Row(Row, Data.Get_Data()) : false );
1191}
1192
1193//---------------------------------------------------------
1194bool 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
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//---------------------------------------------------------
1224bool 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
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//---------------------------------------------------------
1350CSG_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 {
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//---------------------------------------------------------
1404bool 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//---------------------------------------------------------
1423bool 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//---------------------------------------------------------
1443{
1444 return( Create(Matrix) );
1445}
1446
1447//---------------------------------------------------------
1448bool 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//---------------------------------------------------------
1467bool 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//---------------------------------------------------------
1505bool 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//---------------------------------------------------------
1576bool 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
1699CSG_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//---------------------------------------------------------
1756bool 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//---------------------------------------------------------
1845CSG_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//---------------------------------------------------------
1866CSG_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//---------------------------------------------------------
1885CSG_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//---------------------------------------------------------
1910bool 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//---------------------------------------------------------
1933bool 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//---------------------------------------------------------
1948bool 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//---------------------------------------------------------
2053bool 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//---------------------------------------------------------
bool SG_UI_Process_Get_Okay(bool bBlink)
bool SG_UI_Process_Set_Progress(int Position, int Range)
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
signed long long sLong
Definition api_core.h:158
#define SG_T(s)
Definition api_core.h:537
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
#define SG_Char
Definition api_core.h:536
void * Get_Array(void) const
Definition api_core.h:336
sLong Get_NRows(void) const
Definition mat_tools.h:525
bool from_String(const CSG_String &String)
CSG_Matrix Get_Transpose(void) const
double ** Get_Data(void) const
Definition mat_tools.h:528
bool Del_Rows(sLong nRows)
bool Set_Transpose(void)
CSG_Matrix & operator-=(double Scalar)
bool Add_Row(const double *Data=NULL)
int Get_NX(void) const
Definition mat_tools.h:522
bool Ins_Row(sLong Row, const double *Data=NULL)
CSG_Matrix operator*(double Scalar) const
virtual ~CSG_Matrix(void)
bool Destroy(void)
bool Set_Rows(sLong nRows)
bool Create(const CSG_Matrix &Matrix)
bool Set_Size(sLong nRows, sLong nCols)
bool Set_Cols(sLong nCols)
double Get_Determinant(void) const
CSG_Matrix(void)
bool Del_Row(sLong Row)
bool Subtract(const CSG_Matrix &Matrix)
bool Set_Identity(void)
bool Multiply(double Scalar)
bool Ins_Col(sLong Col, const double *Data=NULL)
CSG_Matrix operator-(double Scalar) const
bool operator==(const CSG_Matrix &Matrix) const
CSG_Vector Get_Col(sLong Col) const
int Get_NY(void) const
Definition mat_tools.h:523
CSG_Matrix & operator*=(double Scalar)
bool is_Square(void) const
Definition mat_tools.h:538
CSG_Matrix & operator+=(double Scalar)
bool is_Equal(const CSG_Matrix &Matrix) const
bool Assign(double Scalar)
bool Del_Cols(sLong nCols)
bool Set_Inverse(bool bSilent=true, int nSubSquare=0)
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
bool Add_Rows(sLong nRows)
bool Set_Zero(void)
bool Set_Col(sLong Col, const double *Data)
bool Add_Cols(sLong nCols)
bool Add_Col(const double *Data=NULL)
CSG_Matrix operator+(double Scalar) const
CSG_Matrix Get_Inverse(bool bSilent=true, int nSubSquare=0) const
bool Del_Col(sLong Col)
CSG_Matrix & operator=(double Scalar)
bool Add(double Scalar)
bool Set_Row(sLong Row, const double *Data)
CSG_Vector Get_Row(sLong Row) const
CSG_String Get_Next_Token(void)
bool Has_More_Tokens(void) const
CSG_String AfterFirst(char Character) const
static CSG_String Format(const char *Format,...)
double asDouble(void) const
CSG_Vector(void)
CSG_Vector Get_Scalar_Product(const CSG_Vector &Vector) const
CSG_Vector operator-(double Scalar) const
bool Multiply_Cross(const CSG_Vector &Vector)
bool Add_Rows(sLong nRows)
double Multiply_Scalar(const CSG_Vector &Vector) const
bool Add_Row(double Value=0.)
bool Create(const CSG_Vector &Vector)
bool Del_Row(sLong Row=-1)
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
bool Sort(bool bAscending=true)
bool Set_Rows(sLong nRows)
sLong Get_Size(void) const
Definition mat_tools.h:382
CSG_Vector & operator-=(double Scalar)
CSG_Vector Get_Cross_Product(const CSG_Vector &Vector) const
CSG_Vector & operator*=(double Scalar)
bool is_Equal(const CSG_Vector &Vector) const
bool Flip_Values(void)
bool Multiply(double Scalar)
bool Destroy(void)
int Get_N(void) const
Definition mat_tools.h:384
virtual ~CSG_Vector(void)
bool is_Collinear(const CSG_Vector &Vector) const
double * Get_Data(void) const
Definition mat_tools.h:386
bool Assign(double Scalar)
bool Set_Zero(void)
CSG_Vector operator+(double Scalar) const
double Get_Angle(const CSG_Vector &Vector) const
CSG_Vector & operator+=(double Scalar)
bool Set_Unity(void)
bool from_String(const CSG_String &String)
bool is_Null(void) const
bool Add(double Scalar)
bool Subtract(const CSG_Vector &Vector)
CSG_Vector operator*(double Scalar) const
bool Del_Rows(sLong nRows)
CSG_Vector & operator=(double Scalar)
CSG_Vector Get_Unity(void) const
double Get_Length(void) const
#define B
#define A
bool SG_Matrix_Eigen_Reduction(const CSG_Matrix &Matrix, CSG_Matrix &Eigen_Vectors, CSG_Vector &Eigen_Values, bool bSilent)
bool SG_VectorR3_Rotate(double Vector[3], size_t Axis, double Angle)
CSG_Matrix SG_Matrix_Get_Rotation(double R, bool bDegree)
bool SG_Matrix_LU_Solve(int n, const int *Permutation, const double **Matrix, double *Vector, bool bSilent)
bool SG_Matrix_Tridiagonal_QL(CSG_Matrix &Q, CSG_Vector &d, CSG_Vector &e)
bool SG_Matrix_Triangular_Decomposition(CSG_Matrix &A, CSG_Vector &d, CSG_Vector &e)
bool SG_VectorR2_Rotate(double &x, double &y, double Angle)
bool SG_Matrix_LU_Decomposition(int n, int *Permutation, double **Matrix, bool bSilent, int *nRowChanges)
CSG_Vector operator*(double Scalar, const CSG_Vector &Vector)
bool SG_Matrix_Solve(CSG_Matrix &Matrix, CSG_Vector &Vector, bool bSilent)
#define M_FLT_EPSILON
Definition mat_tools.h:121
SAGA_API_DLL_EXPORT CSG_String SG_Get_Double_asString(double Number, int Width=-1, int Precision=-1, bool bScientific=false)
SAGA_API_DLL_EXPORT int SG_Get_Digit_Count(int Number)
SAGA_API_DLL_EXPORT bool SG_Matrix_LU_Decomposition(int n, int *Permutation, double **Matrix, bool bSilent=true, int *nRowChanges=NULL)
SAGA_API_DLL_EXPORT int SG_Compare_Double(const void *a, const void *b)
#define M_DEG_TO_RAD
Definition mat_tools.h:109
#define M_SET_SIGN(x, sign)
Definition mat_tools.h:145
SAGA_API_DLL_EXPORT bool SG_Matrix_LU_Solve(int n, const int *Permutation, const double **Matrix, double *Vector, bool bSilent=true)