CSL 5.2

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

matrix.hpp
Go to the documentation of this file.
1 ////////////////////////////////
2 // Matrix TCL Lite v1.13
3 // Copyright (c) 1997-2002 Techsoft Pvt. Ltd. (See License.Txt file.)
4 //
5 // Matrix.h: Matrix C++ template class include file
6 // Web: http://www.techsoftpl.com/matrix/
7 // Email: matrix@techsoftpl.com
8 //
9
10 //////////////////////////////
11 // Installation:
12 //
13 // Copy this "matrix.h" file into include directory of your compiler.
14 //
15
16 //////////////////////////////
17 // Note: This matrix template class defines majority of the matrix
18 // operations as overloaded operators or methods. It is assumed that
19 // users of this class is familiar with matrix algebra. We have not
20 // defined any specialization of this template here, so all the instances
21 // of matrix will be created implicitly by the compiler. The data types
22 // tested with this class are float, double, long double, complex<float>,
23 // complex<double> and complex<long double>. Note that this class is not
24 // optimized for performance.
25 //
26 // Since implementation of exception, namespace and template are still
27 // not standardized among the various (mainly old) compilers, you may
28 // encounter compilation error with some compilers. In that case remove
29 // any of the above three features by defining the following macros:
30 //
31 // _NO_NAMESPACE: Define this macro to remove namespace support.
32 //
33 // _NO_EXCEPTION: Define this macro to remove exception handling
34 // and use old style of error handling using function.
35 //
36 // _NO_TEMPLATE: If this macro is defined matrix class of double
37 // type will be generated by default. You can also
38 // generate a different type of matrix like float.
39 //
40 // _SGI_BROKEN_STL: For SGI C++ v.7.2.1 compiler.
41 //
42 // Since all the definitions are also included in this header file as
43 // inline function, some compiler may give warning "inline function
44 // can't be expanded". You may ignore/disable this warning using compiler
45 // switches. All the operators/methods defined in this class have their
46 // natural meaning except the followings:
47 //
48 // Operator/Method Description
49 // --------------- -----------
50 // operator () : This function operator can be used as a
51 // two-dimensional subscript operator to get/set
52 // individual matrix elements.
53 //
54 // operator ! : This operator has been used to calculate inversion
55 // of matrix.
56 //
57 // operator ~ : This operator has been used to return transpose of
58 // a matrix.
59 //
60 // operator ^ : It is used calculate power (by a scalar) of a matrix.
61 // When using this operator in a matrix equation, care
62 // must be taken by parenthesizing it because it has
63 // lower precedence than addition, subtraction,
64 // multiplication and division operators.
65 //
66 // operator >> : It is used to read matrix from input stream as per
67 // standard C++ stream operators.
68 //
69 // operator << : It is used to write matrix to output stream as per
70 // standard C++ stream operators.
71 //
72 // Note that professional version of this package, Matrix TCL Pro 2.11
73 // is optimized for performance and supports many more matrix operations.
74 // It is available from our web site at <http://www.techsoftpl.com/matrix/>.
75 //
76
77 #ifndef __cplusplus
78 #error Must use C++ for the type matrix.
79 #endif
80
81 #if !defined(__STD_MATRIX_H)
82 #define __STD_MATRIX_H
83
84 #define __GXX_WEAK__ 1
85 #define __DEPRECATED 1
86 #define __GNUG__ 4
87 #define __EXCEPTIONS 1
88 #define __private_extern__ extern
89
90 //////////////////////////////
91 // First deal with various shortcomings and incompatibilities of
92 // various (mainly old) versions of popular compilers available.
93 //
94
95 #if defined(__BORLANDC__)
96 #pragma option -w-inl -w-pch
97 #endif
98
99 #if (defined(__BORLANDC__) || _MSC_VER <= 1000 ) && !defined(__GNUG__ )
100 # include <stdio.h>
101 # include <stdlib.h>
102 # include <math.h>
103 # include <iostream.h>
104 # include <string.h>
105 #else
106 # include <string.h>
107 # include <cmath>
108 # include <cstdio>
109 # include <cstdlib>
110 # include <string>
111 # include <iostream>
112 #endif
113
114 #if defined(_MSC_VER) && _MSC_VER <= 1000
115 # define _NO_EXCEPTION // stdexception is not fully supported in MSVC++ 4.0
116 typedef int bool;
117 # if !defined(false)
118 # define false 0
119 # endif
120 # if !defined(true)
121 # define true 1
122 # endif
123 #endif
124
125 #if defined(__BORLANDC__) && !defined(__WIN32__)
126 # define _NO_EXCEPTION // std exception and namespace are not fully
127 # define _NO_NAMESPACE // supported in 16-bit compiler
128 #endif
129
130 #if defined(_MSC_VER) && !defined(_WIN32)
131 # define _NO_EXCEPTION
132 #endif
133
134 #if defined(_NO_EXCEPTION)
135 # define _NO_THROW
136 # define _THROW_MATRIX_ERROR
137 #else
138 # if defined(_MSC_VER)
139 # if _MSC_VER >= 1020
140 # include <stdexcept>
141 # else
142 # include <stdexcpt.h>
143 # endif
144 # elif defined(__MWERKS__)
145 # include <stdexcept>
146 # elif (__GNUC__ >= 2 || (__GNUC__ == 2 && __GNUC_MINOR__ >= 8))
147 # include <stdexcept>
148 # else
149 # include <stdexcep>
150 # endif
151 # define _NO_THROW throw ()
152 # define _THROW_MATRIX_ERROR throw (matrix_error)
153 #endif
154
155 #ifndef __MINMAX_DEFINED
156 # define max(a,b) (((a) > (b)) ? (a) : (b))
157 # define min(a,b) (((a) < (b)) ? (a) : (b))
158 #endif
159
160 #if defined(_MSC_VER)
161 #undef _MSC_EXTENSIONS // To include overloaded abs function definitions!
162 #endif
163
164 #if (defined(__BORLANDC__) || _MSC_VER ) && !defined(__GNUG__ )
165 inline float abs (float v) { return (float)fabs(v); }
166 inline double abs (double v) { return fabs(v); }
167 inline long double abs (long double v) { return fabsl(v); }
168 #endif
169
170 #if defined(__GNUG__) || defined(__MWERKS__) || (defined(__BORLANDC__) && (__BORLANDC__ >= 0x540))
171 #define FRIEND_FUN_TEMPLATE <>
172 #else
173 #define FRIEND_FUN_TEMPLATE
174 #endif
175
176 #if defined(_MSC_VER) && _MSC_VER <= 1020 // MSVC++ 4.0/4.2 does not
177 # define _NO_NAMESPACE // support "std" namespace
178 #endif
179
180 #if !defined(_NO_NAMESPACE)
181 #if defined(_SGI_BROKEN_STL ) // For SGI C++ v.7.2.1 compiler
182 namespace std { }
183 #endif
184 using namespace std;
185 #endif
186
187 #ifndef _NO_NAMESPACE
188 namespace math {
189 #endif
190
191 #if !defined(_NO_EXCEPTION)
192 class matrix_error : public logic_error
193 {
194 public:
195 matrix_error (const string& what_arg) : logic_error(what_arg) { }
196 };
197 #define REPORT_ERROR(ErrormMsg) throw matrix_error(ErrormMsg);
198 #else
199 inline void _matrix_error (const char* pErrMsg)
200 {
201 cerr << pErrMsg << endl;
202 exit(1);
203 }
204 #define REPORT_ERROR(ErrormMsg) _matrix_error(ErrormMsg);
205 #endif
206
207 #if !defined(_NO_TEMPLATE)
208 # define MAT_TEMPLATE template <class T>
209 # define matrixT matrix<T>
210 #else
211 # define MAT_TEMPLATE
212 # define matrixT matrix
213 # ifdef MATRIX_TYPE
214 typedef MATRIX_TYPE T;
215 # else
216 typedef double T;
217 # endif
218 #endif
219
220
221 MAT_TEMPLATE
222 class matrix
223 {
224 public:
225 // Constructors
226 matrix (const matrixT& m);
227 matrix (size_t row = 6, size_t col = 6);
228
229 // Destructor
230 ~matrix ();
231
232 // Assignment operators
233 matrixT& operator = (const matrixT& m) _NO_THROW;
234
235 // Value extraction method
236 size_t RowNo () const { return _m->Row; }
237 size_t ColNo () const { return _m->Col; }
238
239 // Subscript operator
240 T& operator () (size_t row, size_t col) _THROW_MATRIX_ERROR;
241 T operator () (size_t row, size_t col) const _THROW_MATRIX_ERROR;
242
243 // Unary operators
244 matrixT operator + () _NO_THROW { return *this; }
245 matrixT operator - () _NO_THROW;
246
247 // Combined assignment - calculation operators
248 matrixT& operator += (const matrixT& m) _THROW_MATRIX_ERROR;
249 matrixT& operator -= (const matrixT& m) _THROW_MATRIX_ERROR;
250 matrixT& operator *= (const matrixT& m) _THROW_MATRIX_ERROR;
251 matrixT& operator *= (const T& c) _NO_THROW;
252 matrixT& operator /= (const T& c) _NO_THROW;
253 matrixT& operator ^= (const size_t& pow) _THROW_MATRIX_ERROR;
254
255 // Miscellaneous -methods
256 void Null (const size_t& row, const size_t& col) _NO_THROW;
257 void Null () _NO_THROW;
258 void Unit (const size_t& row) _NO_THROW;
259 void Unit () _NO_THROW;
260 void SetSize (size_t row, size_t col) _NO_THROW;
261
262 // Utility methods
263 matrixT Solve (const matrixT& v) const _THROW_MATRIX_ERROR;
264 matrixT Adj () _THROW_MATRIX_ERROR;
265 matrixT Inv () _THROW_MATRIX_ERROR;
266 T Det () const _THROW_MATRIX_ERROR;
267 T Norm () _NO_THROW;
268 T Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR;
269 T Cond () _NO_THROW;
270
271 // Type of matrices
272 bool IsSquare () _NO_THROW { return (_m->Row == _m->Col); }
273 bool IsSingular () _NO_THROW;
274 bool IsDiagonal () _NO_THROW;
275 bool IsScalar () _NO_THROW;
276 bool IsUnit () _NO_THROW;
277 bool IsNull () _NO_THROW;
278 bool IsSymmetric () _NO_THROW;
279 bool IsSkewSymmetric () _NO_THROW;
280 bool IsUpperTriangular () _NO_THROW;
281 bool IsLowerTriangular () _NO_THROW;
282
283 private:
284 struct base_mat
285 {
286 T **Val;
287 size_t Row, Col, RowSiz, ColSiz;
288 int Refcnt;
289
290 base_mat (size_t row, size_t col, T** v)
291 {
292 Row = row; RowSiz = row;
293 Col = col; ColSiz = col;
294 Refcnt = 1;
295
296 Val = new T* [row];
297 size_t rowlen = col * sizeof(T);
298
299 for (size_t i=0; i < row; i++)
300 {
301 Val[i] = new T [col];
302 if (v) memcpy(Val[i], v[i], rowlen);
303 }
304 }
305 ~base_mat ()
306 {
307 for (size_t i=0; i < RowSiz; i++)
308 delete [] Val[i];
309 delete [] Val;
310 }
311 };
312 base_mat *_m;
313
314 void clone ();
315 void realloc (size_t row, size_t col);
316 int pivot (size_t row);
317 };
318
319 #if defined(_MSC_VER) && _MSC_VER <= 1020
320 # undef _NO_THROW // MSVC++ 4.0/4.2 does not support
321 # undef _THROW_MATRIX_ERROR // exception specification in definition
322 # define _NO_THROW
323 # define _THROW_MATRIX_ERROR
324 #endif
325
326 // constructor
327 MAT_TEMPLATE inline
328 matrixT::matrix (size_t row, size_t col)
329 {
330 _m = new base_mat(row, col, 0);
331 }
332
333 // copy constructor
334 MAT_TEMPLATE inline
335 matrixT::matrix (const matrixT& m)
336 {
337 _m = m._m;
338 _m->Refcnt++;
339 }
340
341 // Internal copy constructor
342 MAT_TEMPLATE inline void
343 matrixT::clone ()
344 {
345 _m->Refcnt--;
346 _m = new base_mat(_m->Row, _m->Col, _m->Val);
347 }
348
349 // destructor
350 MAT_TEMPLATE inline
351 matrixT::~matrix ()
352 {
353 if (--_m->Refcnt == 0) delete _m;
354 }
355
356 // assignment operator
357 MAT_TEMPLATE inline matrixT&
358 matrixT::operator = (const matrixT& m) _NO_THROW
359 {
360 m._m->Refcnt++;
361 if (--_m->Refcnt == 0) delete _m;
362 _m = m._m;
363 return *this;
364 }
365
366 // reallocation method
367 MAT_TEMPLATE inline void
368 matrixT::realloc (size_t row, size_t col)
369 {
370 if (row == _m->RowSiz && col == _m->ColSiz)
371 {
372 _m->Row = _m->RowSiz;
373 _m->Col = _m->ColSiz;
374 return;
375 }
376
377 base_mat *m1 = new base_mat(row, col, NULL);
378 size_t colSize = min(_m->Col,col) * sizeof(T);
379 size_t minRow = min(_m->Row,row);
380
381 for (size_t i=0; i < minRow; i++)
382 memcpy(m1->Val[i], _m->Val[i], colSize);
383
384 if (--_m->Refcnt == 0)
385 delete _m;
386 _m = m1;
387
388 return;
389 }
390
391 // public method for resizing matrix
392 MAT_TEMPLATE inline void
393 matrixT::SetSize (size_t row, size_t col) _NO_THROW
394 {
395 size_t i,j;
396 size_t oldRow = _m->Row;
397 size_t oldCol = _m->Col;
398
399 if (row != _m->RowSiz || col != _m->ColSiz)
400 realloc(row, col);
401
402 for (i=oldRow; i < row; i++)
403 for (j=0; j < col; j++)
404 _m->Val[i][j] = T(0);
405
406 for (i=0; i < row; i++)
407 for (j=oldCol; j < col; j++)
408 _m->Val[i][j] = T(0);
409
410 return;
411 }
412
413 // subscript operator to get/set individual elements
414 MAT_TEMPLATE inline T&
415 matrixT::operator () (size_t row, size_t col) _THROW_MATRIX_ERROR
416 {
417 if (row >= _m->Row || col >= _m->Col)
418 REPORT_ERROR("matrixT::operator(): Index out of range!");
419 if (_m->Refcnt > 1) clone();
420 return _m->Val[row][col];
421 }
422
423 // subscript operator to get/set individual elements
424 MAT_TEMPLATE inline T
425 matrixT::operator () (size_t row, size_t col) const _THROW_MATRIX_ERROR
426 {
427 if (row >= _m->Row || col >= _m->Col)
428 REPORT_ERROR("matrixT::operator(): Index out of range!");
429 return _m->Val[row][col];
430 }
431
432 // input stream function
433 MAT_TEMPLATE inline istream&
434 operator >> (istream& istrm, matrixT& m)
435 {
436 for (size_t i=0; i < m.RowNo(); i++)
437 for (size_t j=0; j < m.ColNo(); j++)
438 {
439 T x;
440 istrm >> x;
441 m(i,j) = x;
442 }
443 return istrm;
444 }
445
446 // output stream function
447 MAT_TEMPLATE inline ostream&
448 operator << (ostream& ostrm, const matrixT& m)
449 {
450 for (size_t i=0; i < m.RowNo(); i++)
451 {
452 for (size_t j=0; j < m.ColNo(); j++)
453 {
454 T x = m(i,j);
455 ostrm << x << '\t';
456 }
457 ostrm << endl;
458 }
459 return ostrm;
460 }
461
462
463 // logical equal-to operator
464 MAT_TEMPLATE inline bool
465 operator == (const matrixT& m1, const matrixT& m2) _NO_THROW
466 {
467 if (m1.RowNo() != m2.RowNo() || m1.ColNo() != m2.ColNo())
468 return false;
469
470 for (size_t i=0; i < m1.RowNo(); i++)
471 for (size_t j=0; j < m1.ColNo(); j++)
472 if (m1(i,j) != m2(i,j))
473 return false;
474
475 return true;
476 }
477
478 // logical no-equal-to operator
479 MAT_TEMPLATE inline bool
480 operator != (const matrixT& m1, const matrixT& m2) _NO_THROW
481 {
482 return (m1 == m2) ? false : true;
483 }
484
485 // combined addition and assignment operator
486 MAT_TEMPLATE inline matrixT&
487 matrixT::operator += (const matrixT& m) _THROW_MATRIX_ERROR
488 {
489 if (_m->Row != m._m->Row || _m->Col != m._m->Col)
490 REPORT_ERROR("matrixT::operator+= : Inconsistent matrix sizes in addition!");
491 if (_m->Refcnt > 1) clone();
492 for (size_t i=0; i < m._m->Row; i++)
493 for (size_t j=0; j < m._m->Col; j++)
494 _m->Val[i][j] += m._m->Val[i][j];
495 return *this;
496 }
497
498 // combined subtraction and assignment operator
499 MAT_TEMPLATE inline matrixT&
500 matrixT::operator -= (const matrixT& m) _THROW_MATRIX_ERROR
501 {
502 if (_m->Row != m._m->Row || _m->Col != m._m->Col)
503 REPORT_ERROR("matrixT::operator-= : Inconsistent matrix sizes in subtraction!");
504 if (_m->Refcnt > 1) clone();
505 for (size_t i=0; i < m._m->Row; i++)
506 for (size_t j=0; j < m._m->Col; j++)
507 _m->Val[i][j] -= m._m->Val[i][j];
508 return *this;
509 }
510
511 // combined scalar multiplication and assignment operator
512 MAT_TEMPLATE inline matrixT&
513 matrixT::operator *= (const T& c) _NO_THROW
514 {
515 if (_m->Refcnt > 1) clone();
516 for (size_t i=0; i < _m->Row; i++)
517 for (size_t j=0; j < _m->Col; j++)
518 _m->Val[i][j] *= c;
519 return *this;
520 }
521
522 // combined matrix multiplication and assignment operator
523 MAT_TEMPLATE inline matrixT&
524 matrixT::operator *= (const matrixT& m) _THROW_MATRIX_ERROR
525 {
526 if (_m->Col != m._m->Row)
527 REPORT_ERROR("matrixT::operator*= : Inconsistent matrix sizes in multiplication!");
528
529 matrixT temp(_m->Row,m._m->Col);
530
531 for (size_t i=0; i < _m->Row; i++)
532 for (size_t j=0; j < m._m->Col; j++)
533 {
534 temp._m->Val[i][j] = T(0);
535 for (size_t k=0; k < _m->Col; k++)
536 temp._m->Val[i][j] += _m->Val[i][k] * m._m->Val[k][j];
537 }
538 *this = temp;
539
540 return *this;
541 }
542
543 // combined scalar division and assignment operator
544 MAT_TEMPLATE inline matrixT&
545 matrixT::operator /= (const T& c) _NO_THROW
546 {
547 if (_m->Refcnt > 1) clone();
548 for (size_t i=0; i < _m->Row; i++)
549 for (size_t j=0; j < _m->Col; j++)
550 _m->Val[i][j] /= c;
551
552 return *this;
553 }
554
555 // combined power and assignment operator
556 MAT_TEMPLATE inline matrixT&
557 matrixT::operator ^= (const size_t& pow) _THROW_MATRIX_ERROR
558 {
559 matrixT temp(*this);
560
561 for (size_t i=2; i <= pow; i++)
562 *this = *this * temp;
563
564 return *this;
565 }
566
567 // unary negation operator
568 MAT_TEMPLATE inline matrixT
569 matrixT::operator - () _NO_THROW
570 {
571 matrixT temp(_m->Row,_m->Col);
572
573 for (size_t i=0; i < _m->Row; i++)
574 for (size_t j=0; j < _m->Col; j++)
575 temp._m->Val[i][j] = - _m->Val[i][j];
576
577 return temp;
578 }
579
580 // binary addition operator
581 MAT_TEMPLATE inline matrixT
582 operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
583 {
584 matrixT temp = m1;
585 temp += m2;
586 return temp;
587 }
588
589 // binary subtraction operator
590 MAT_TEMPLATE inline matrixT
591 operator - (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
592 {
593 matrixT temp = m1;
594 temp -= m2;
595 return temp;
596 }
597
598 // binary scalar multiplication operator
599 MAT_TEMPLATE inline matrixT
600 operator * (const matrixT& m, const T& no) _NO_THROW
601 {
602 matrixT temp = m;
603 temp *= no;
604 return temp;
605 }
606
607
608 // binary scalar multiplication operator
609 MAT_TEMPLATE inline matrixT
610 operator * (const T& no, const matrixT& m) _NO_THROW
611 {
612 return (m * no);
613 }
614
615 // binary matrix element multiplication operator
616 MAT_TEMPLATE inline matrixT
617 operator * (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
618 {
619 matrixT temp = m1;
620 temp *= m2;
621 return temp;
622 }
623
624
625 // binary scalar division operator
626 MAT_TEMPLATE inline matrixT
627 operator / (const matrixT& m, const T& no) _NO_THROW
628 {
629 return (m * (T(1) / no));
630 }
631
632
633 // binary scalar division operator
634 MAT_TEMPLATE inline matrixT
635 operator / (const T& no, const matrixT& m) _THROW_MATRIX_ERROR
636 {
637 return ( ! m * no);
638 }
639
640 // binary matrix division operator
641 MAT_TEMPLATE inline matrixT
642 operator / (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
643 {
644 return (m1 * !m2);
645 }
646
647 // binary power operator
648 MAT_TEMPLATE inline matrixT
649 operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR
650 {
651 matrixT temp = m;
652 temp ^= pow;
653 return temp;
654 }
655
656 // unary transpose operator
657 MAT_TEMPLATE inline matrixT
658 operator ~ (const matrixT& m) _NO_THROW
659 {
660 matrixT temp(m.ColNo(),m.RowNo());
661
662 for (size_t i=0; i < m.RowNo(); i++)
663 for (size_t j=0; j < m.ColNo(); j++)
664 {
665 T x = m(i,j);
666 temp(j,i) = x;
667 }
668 return temp;
669 }
670
671 // unary inversion operator
672 MAT_TEMPLATE inline matrixT
673 operator ! (const matrixT m) _THROW_MATRIX_ERROR
674 {
675 matrixT temp = m;
676 return temp.Inv();
677 }
678
679 // inversion function
680 MAT_TEMPLATE inline matrixT
681 matrixT::Inv () _THROW_MATRIX_ERROR
682 {
683 size_t i,j,k;
684 T a1,a2,*rowptr;
685
686 if (_m->Row != _m->Col)
687 REPORT_ERROR("matrixT::operator!: Inversion of a non-square matrix");
688
689 matrixT temp(_m->Row,_m->Col);
690 if (_m->Refcnt > 1) clone();
691
692
693 temp.Unit();
694 for (k=0; k < _m->Row; k++)
695 {
696 int indx = pivot(k);
697 if (indx == -1)
698 REPORT_ERROR("matrixT::operator!: Inversion of a singular matrix");
699
700 if (indx != 0)
701 {
702 rowptr = temp._m->Val[k];
703 temp._m->Val[k] = temp._m->Val[indx];
704 temp._m->Val[indx] = rowptr;
705 }
706 a1 = _m->Val[k][k];
707 for (j=0; j < _m->Row; j++)
708 {
709 _m->Val[k][j] /= a1;
710 temp._m->Val[k][j] /= a1;
711 }
712 for (i=0; i < _m->Row; i++)
713 if (i != k)
714 {
715 a2 = _m->Val[i][k];
716 for (j=0; j < _m->Row; j++)
717 {
718 _m->Val[i][j] -= a2 * _m->Val[k][j];
719 temp._m->Val[i][j] -= a2 * temp._m->Val[k][j];
720 }
721 }
722 }
723 return temp;
724 }
725
726 // solve simultaneous equation
727 MAT_TEMPLATE inline matrixT
728 matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR
729 {
730 size_t i,j,k;
731 T a1;
732
733 if ( ! (_m->Row == _m->Col && _m->Col == v._m->Row))
734 REPORT_ERROR("matrixT::Solve():Inconsistent matrices!");
735
736 matrixT temp(_m->Row,_m->Col+v._m->Col);
737 for (i=0; i < _m->Row; i++)
738 {
739 for (j=0; j < _m->Col; j++)
740 temp._m->Val[i][j] = _m->Val[i][j];
741 for (k=0; k < v._m->Col; k++)
742 temp._m->Val[i][_m->Col+k] = v._m->Val[i][k];
743 }
744 for (k=0; k < _m->Row; k++)
745 {
746 int indx = temp.pivot(k);
747 if (indx == -1)
748 REPORT_ERROR("matrixT::Solve(): Singular matrix!");
749
750 a1 = temp._m->Val[k][k];
751 for (j=k; j < temp._m->Col; j++)
752 temp._m->Val[k][j] /= a1;
753
754 for (i=k+1; i < _m->Row; i++)
755 {
756 a1 = temp._m->Val[i][k];
757 for (j=k; j < temp._m->Col; j++)
758 temp._m->Val[i][j] -= a1 * temp._m->Val[k][j];
759 }
760 }
761 matrixT s(v._m->Row,v._m->Col);
762 for (k=0; k < v._m->Col; k++)
763 for (int m=int(_m->Row)-1; m >= 0; m--)
764 {
765 s._m->Val[m][k] = temp._m->Val[m][_m->Col+k];
766 for (j=m+1; j < _m->Col; j++)
767 s._m->Val[m][k] -= temp._m->Val[m][j] * s._m->Val[j][k];
768 }
769 return s;
770 }
771
772 // set zero to all elements of this matrix
773 MAT_TEMPLATE inline void
774 matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
775 {
776 if (row != _m->Row || col != _m->Col)
777 realloc(row,col);
778
779 if (_m->Refcnt > 1)
780 clone();
781
782 for (size_t i=0; i < _m->Row; i++)
783 for (size_t j=0; j < _m->Col; j++)
784 _m->Val[i][j] = T(0);
785 return;
786 }
787
788 // set zero to all elements of this matrix
789 MAT_TEMPLATE inline void
790 matrixT::Null() _NO_THROW
791 {
792 if (_m->Refcnt > 1) clone();
793 for (size_t i=0; i < _m->Row; i++)
794 for (size_t j=0; j < _m->Col; j++)
795 _m->Val[i][j] = T(0);
796 return;
797 }
798
799 // set this matrix to unity
800 MAT_TEMPLATE inline void
801 matrixT::Unit (const size_t& row) _NO_THROW
802 {
803 if (row != _m->Row || row != _m->Col)
804 realloc(row, row);
805
806 if (_m->Refcnt > 1)
807 clone();
808
809 for (size_t i=0; i < _m->Row; i++)
810 for (size_t j=0; j < _m->Col; j++)
811 _m->Val[i][j] = i == j ? T(1) : T(0);
812 return;
813 }
814
815 // set this matrix to unity
816 MAT_TEMPLATE inline void
817 matrixT::Unit () _NO_THROW
818 {
819 if (_m->Refcnt > 1) clone();
820 size_t row = min(_m->Row,_m->Col);
821 _m->Row = _m->Col = row;
822
823 for (size_t i=0; i < _m->Row; i++)
824 for (size_t j=0; j < _m->Col; j++)
825 _m->Val[i][j] = i == j ? T(1) : T(0);
826 return;
827 }
828
829 // private partial pivoting method
830 MAT_TEMPLATE inline int
831 matrixT::pivot (size_t row)
832 {
833 int k = int(row);
834 double amax,temp;
835
836 amax = -1;
837 for (size_t i=row; i < _m->Row; i++)
838 if ((temp = abs(_m->Val[i][row])) > amax && temp != 0.0)
839 {
840 amax = temp;
841 k = i;
842 }
843 if (_m->Val[k][row] == T(0))
844 return -1;
845 if (k != int(row))
846 {
847 T* rowptr = _m->Val[k];
848 _m->Val[k] = _m->Val[row];
849 _m->Val[row] = rowptr;
850 return k;
851 }
852 return 0;
853 }
854
855 // calculate the determinant of a matrix
856 MAT_TEMPLATE inline T
857 matrixT::Det () const _THROW_MATRIX_ERROR
858 {
859 size_t i,j,k;
860 T piv,detVal = T(1);
861
862 if (_m->Row != _m->Col)
863 REPORT_ERROR("matrixT::Det(): Determinant a non-square matrix!");
864
865 matrixT temp(*this);
866 if (temp._m->Refcnt > 1) temp.clone();
867
868 for (k=0; k < _m->Row; k++)
869 {
870 int indx = temp.pivot(k);
871 if (indx == -1)
872 return 0;
873 if (indx != 0)
874 detVal = - detVal;
875 detVal = detVal * temp._m->Val[k][k];
876 for (i=k+1; i < _m->Row; i++)
877 {
878 piv = temp._m->Val[i][k] / temp._m->Val[k][k];
879 for (j=k+1; j < _m->Row; j++)
880 temp._m->Val[i][j] -= piv * temp._m->Val[k][j];
881 }
882 }
883 return detVal;
884 }
885
886 // calculate the norm of a matrix
887 MAT_TEMPLATE inline T
888 matrixT::Norm () _NO_THROW
889 {
890 T retVal = T(0);
891
892 for (size_t i=0; i < _m->Row; i++)
893 for (size_t j=0; j < _m->Col; j++)
894 retVal += _m->Val[i][j] * _m->Val[i][j];
895 retVal = sqrt(retVal);
896
897 return retVal;
898 }
899
900 // calculate the condition number of a matrix
901 MAT_TEMPLATE inline T
902 matrixT::Cond () _NO_THROW
903 {
904 matrixT inv = ! (*this);
905 return (Norm() * inv.Norm());
906 }
907
908 // calculate the cofactor of a matrix for a given element
909 MAT_TEMPLATE inline T
910 matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR
911 {
912 size_t i,i1,j,j1;
913
914 if (_m->Row != _m->Col)
915 REPORT_ERROR("matrixT::Cofact(): Cofactor of a non-square matrix!");
916
917 if (row > _m->Row || col > _m->Col)
918 REPORT_ERROR("matrixT::Cofact(): Index out of range!");
919
920 matrixT temp (_m->Row-1,_m->Col-1);
921
922 for (i=i1=0; i < _m->Row; i++)
923 {
924 if (i == row)
925 continue;
926 for (j=j1=0; j < _m->Col; j++)
927 {
928 if (j == col)
929 continue;
930 temp._m->Val[i1][j1] = _m->Val[i][j];
931 j1++;
932 }
933 i1++;
934 }
935 T cof = temp.Det();
936 if ((row+col)%2 == 1)
937 cof = -cof;
938
939 return cof;
940 }
941
942
943 // calculate adjoin of a matrix
944 MAT_TEMPLATE inline matrixT
945 matrixT::Adj () _THROW_MATRIX_ERROR
946 {
947 if (_m->Row != _m->Col)
948 REPORT_ERROR("matrixT::Adj(): Adjoin of a non-square matrix.");
949
950 matrixT temp(_m->Row,_m->Col);
951
952 for (size_t i=0; i < _m->Row; i++)
953 for (size_t j=0; j < _m->Col; j++)
954 temp._m->Val[j][i] = Cofact(i,j);
955 return temp;
956 }
957
958 // Determine if the matrix is singular
959 MAT_TEMPLATE inline bool
960 matrixT::IsSingular () _NO_THROW
961 {
962 if (_m->Row != _m->Col)
963 return false;
964 return (Det() == T(0));
965 }
966
967 // Determine if the matrix is diagonal
968 MAT_TEMPLATE inline bool
969 matrixT::IsDiagonal () _NO_THROW
970 {
971 if (_m->Row != _m->Col)
972 return false;
973 for (size_t i=0; i < _m->Row; i++)
974 for (size_t j=0; j < _m->Col; j++)
975 if (i != j && _m->Val[i][j] != T(0))
976 return false;
977 return true;
978 }
979
980 // Determine if the matrix is scalar
981 MAT_TEMPLATE inline bool
982 matrixT::IsScalar () _NO_THROW
983 {
984 if ( ! IsDiagonal())
985 return false;
986 T v = _m->Val[0][0];
987 for (size_t i=1; i < _m->Row; i++)
988 if (_m->Val[i][i] != v)
989 return false;
990 return true;
991 }
992
993 // Determine if the matrix is a unit matrix
994 MAT_TEMPLATE inline bool
995 matrixT::IsUnit () _NO_THROW
996 {
997 if (IsScalar() && _m->Val[0][0] == T(1))
998 return true;
999 return false;
1000 }
1001
1002 // Determine if this is a null matrix
1003 MAT_TEMPLATE inline bool
1004 matrixT::IsNull () _NO_THROW
1005 {
1006 for (size_t i=0; i < _m->Row; i++)
1007 for (size_t j=0; j < _m->Col; j++)
1008 if (_m->Val[i][j] != T(0))
1009 return false;
1010 return true;
1011 }
1012
1013 // Determine if the matrix is symmetric
1014 MAT_TEMPLATE inline bool
1015 matrixT::IsSymmetric () _NO_THROW
1016 {
1017 if (_m->Row != _m->Col)
1018 return false;
1019 for (size_t i=0; i < _m->Row; i++)
1020 for (size_t j=0; j < _m->Col; j++)
1021 if (_m->Val[i][j] != _m->Val[j][i])
1022 return false;
1023 return true;
1024 }
1025
1026 // Determine if the matrix is skew-symmetric
1027 MAT_TEMPLATE inline bool
1028 matrixT::IsSkewSymmetric () _NO_THROW
1029 {
1030 if (_m->Row != _m->Col)
1031 return false;
1032 for (size_t i=0; i < _m->Row; i++)
1033 for (size_t j=0; j < _m->Col; j++)
1034 if (_m->Val[i][j] != -_m->Val[j][i])
1035 return false;
1036 return true;
1037 }
1038
1039 // Determine if the matrix is upper triangular
1040 MAT_TEMPLATE inline bool
1041 matrixT::IsUpperTriangular () _NO_THROW
1042 {
1043 if (_m->Row != _m->Col)
1044 return false;
1045 for (size_t i=1; i < _m->Row; i++)
1046 for (size_t j=0; j < i-1; j++)
1047 if (_m->Val[i][j] != T(0))
1048 return false;
1049 return true;
1050 }
1051
1052 // Determine if the matrix is lower triangular
1053 MAT_TEMPLATE inline bool
1054 matrixT::IsLowerTriangular () _NO_THROW
1055 {
1056 if (_m->Row != _m->Col)
1057 return false;
1058
1059 for (size_t j=1; j < _m->Col; j++)
1060 for (size_t i=0; i < j-1; i++)
1061 if (_m->Val[i][j] != T(0))
1062 return false;
1063
1064 return true;
1065 }
1066
1067 #ifndef _NO_NAMESPACE
1068 }
1069 #endif
1070
1071 #endif //__STD_MATRIX_H