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