- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

matrix.h
説明を見る。
1 /****************************************************************************************************************************
2 **
3 ** MIST ( Media Integration Standard Toolkit )
4 **
5 ** matrix template class implementation using expression template technique.
6 **
7 **
8 ** All matrix operations are described in row order.
9 ** All matrix elements are allocated as one dimensional array on the memory space.
10 **
11 ** ex) The order of 3x3 matrix are
12 **
13 ** 1 4 7
14 ** 2 5 8
15 ** 3 6 9
16 **
17 **
18 ** Please use matrix( row, col ) accesses operation.
19 **
20 **
21 **
22 ** We developed these programs since 2003/09/05.
23 **
24 ** $LastChangedDate:: 2011-04-26 13:10:24 #$
25 ** $LastChangedRevision: 1365 $
26 ** $LastChangedBy: ddeguchi $
27 ** $HeadURL: http://localhost/svn/mist/trunk/mist/matrix.h $
28 **
29 ** Copyright ***********************.
30 ** All Rights Reserved.
31 **
32 ****************************************************************************************************************************/
33
34 //
35 // Copyright (c) 2003-2011, MIST Project, Nagoya University
36 // All rights reserved.
37 //
38 // Redistribution and use in source and binary forms, with or without modification,
39 // are permitted provided that the following conditions are met:
40 //
41 // 1. Redistributions of source code must retain the above copyright notice,
42 // this list of conditions and the following disclaimer.
43 //
44 // 2. Redistributions in binary form must reproduce the above copyright notice,
45 // this list of conditions and the following disclaimer in the documentation
46 // and/or other materials provided with the distribution.
47 //
48 // 3. Neither the name of the Nagoya University nor the names of its contributors
49 // may be used to endorse or promote products derived from this software
50 // without specific prior written permission.
51 //
52 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
53 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
54 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
55 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
56 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
57 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
58 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
59 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
60 //
61
62
63
68 #ifndef __INCLUDE_MIST_MATRIX__
69 #define __INCLUDE_MIST_MATRIX__
70
71 #include <complex>
72
73 #ifndef __INCLUDE_MIST_H__
74 #include "mist.h"
75 #endif
76
77 #ifdef _OPENMP
78 #include <omp.h>
79 #endif
80
81 #include <cmath>
82
83
84 // mist名前空間の始まり
85 _MIST_BEGIN
86
87
88
89 namespace __numeric__
90 {
91 template < class T >
92 struct is_complex
93 {
94 _MIST_CONST( bool, value, false );
95 };
96
97 #if defined(__MIST_MSVC__) && __MIST_MSVC__ < 7
98
99 #define IS_COMPLEX( type ) \
100 template < >\
101 struct is_complex< std::complex< type > >\
102 {\
103 enum{ value = true };\
104 };\
105
106 // 各型に対する特殊化
107 IS_COMPLEX(unsigned char)
108 IS_COMPLEX(unsigned short)
109 IS_COMPLEX(unsigned int)
110 IS_COMPLEX(unsigned long)
111 IS_COMPLEX(signed char)
112 IS_COMPLEX(signed short)
113 IS_COMPLEX(signed int)
114 IS_COMPLEX(signed long)
115 IS_COMPLEX(bool)
116 IS_COMPLEX(char)
117 IS_COMPLEX(float)
118 IS_COMPLEX(double)
119 IS_COMPLEX(long double)
120
121 #undef IS_COLOR
122
123 #else
124
125 template < class T >
126 struct is_complex< std::complex< T > >
127 {
128 _MIST_CONST( bool, value, true );
129 };
130
131 #endif
132
133 template< bool b >
134 struct value_compare
135 {
136 template < class T > static bool le( const T &v1, const T &v2 ){ return( v1 <= v2 ); }
137 template < class T > static bool ge( const T &v1, const T &v2 ){ return( v1 >= v2 ); }
138 template < class T > static bool lt( const T &v1, const T &v2 ){ return( v1 < v2 ); }
139 template < class T > static bool gt( const T &v1, const T &v2 ){ return( v1 > v2 ); }
140 template < class T > static bool eq( const T &v1, const T &v2 ){ return( v1 == v2 ); }
141 template < class T > static bool eq( const T &v1, const T &v2, const double delta ){ return( std::abs( v1 - v2 ) < delta ); }
142 template < class T > static bool is_zero( const T &v ){ return( v == T( ) ); }
143 };
144
145 template< >
146 struct value_compare< true >
147 {
148 template < class T >
149 static bool lt( const T &v1, const T &v2 )
150 {
151 if( v1.real( ) < v2.real( ) )
152 {
153 return( true );
154 }
155 else if( v1.real( ) > v2.real( ) )
156 {
157 return( false );
158 }
159 else if( v1.imag( ) < v2.imag( ) )
160 {
161 return( true );
162 }
163 else
164 {
165 return( false );
166 }
167 }
168 template < class T > static bool ge( const T &v1, const T &v2 ){ return( !lt( v1, v2 ) ); }
169 template < class T > static bool le( const T &v1, const T &v2 ){ return( !lt( v2, v1 ) ); }
170 template < class T > static bool gt( const T &v1, const T &v2 ){ return( lt( v2, v1 ) ); }
171 template < class T > static bool eq( const T &v1, const T &v2 )
172 {
173 return( v1.real( ) == v2.real( ) && v1.imag( ) == v2.imag( ) );
174 }
175 template < class T > static bool eq( const T &v1, const T &v2, const double delta )
176 {
177 return( std::abs( v1.real( ) - v2.real( ) ) < delta && std::abs( v1.imag( ) - v2.imag( ) ) < delta );
178 }
179 template < class T > static bool is_zero( const T &v ){ return( v.real( ) == 0 && v.imag( ) == 0 ); }
180 };
181
182 // 共役複素数を返す関数
183 template< class T > std::complex< T > conjugate( const std::complex< T > &c ) { return std::conj( c ); }
184 template< class T > const T &conjugate( const T &c ) { return c; }
185 }
186
187
196 template < class T, class Allocator = ::std::allocator< T > >
197 class matrix : public array< T, Allocator >
198 {
199 public:
200 typedef Allocator allocator_type;
201 typedef typename Allocator::reference reference;
202 typedef typename Allocator::const_reference const_reference;
203 typedef typename Allocator::value_type value_type;
204 typedef typename Allocator::size_type size_type;
205 typedef typename Allocator::difference_type difference_type;
206 typedef typename Allocator::pointer pointer;
207 typedef typename Allocator::const_pointer const_pointer;
208
210 typedef mist_iterator1< T, ptrdiff_t, pointer, reference > iterator;
211
213 typedef mist_iterator1< T, ptrdiff_t, const_pointer, const_reference > const_iterator;
214
216 typedef mist_reverse_iterator< mist_iterator1< T, ptrdiff_t, pointer, reference > > reverse_iterator;
217
219 typedef mist_reverse_iterator< mist_iterator1< T, ptrdiff_t, const_pointer, const_reference > > const_reverse_iterator;
220
222 template < class TT, class AAllocator = std::allocator< TT > >
223 struct rebind
224 {
225 typedef matrix< TT, AAllocator > other;
226 };
227
228 private:
229 typedef array< T, Allocator > base;
230 size_type size2_;
231 size_type size1_;
232
233 private:
234 // rows や cols との混乱を避けるため，matrix では使えないようにする
235 size_type width( ) const;
236 size_type height( ) const;
237 size_type depth( ) const;
238
239 public:
251 bool resize( size_type num )
252 {
253 if( base::resize( num ) )
254 {
255 size1_ = num;
256 size2_ = 1;
257 return( true );
258 }
259 else
260 {
261 size1_ = size2_ = 0;
262 return( false );
263 }
264 }
265
266
279 bool resize( size_type nrows, size_type ncols )
280 {
281 if( base::resize( nrows * ncols ) )
282 {
283 size1_ = nrows;
284 size2_ = ncols;
285 return( true );
286 }
287 else
288 {
289 size1_ = size2_ = 0;
290 return( false );
291 }
292 }
293
305 bool resize( size_type nrows, size_type ncols, const T &val )
306 {
307 if( base::resize( nrows * ncols, val ) )
308 {
309 size1_ = nrows;
310 size2_ = ncols;
311 return( true );
312 }
313 else
314 {
315 size1_ = size2_ = 0;
316 return( false );
317 }
318 }
319
331 bool trim( matrix &out, size_type row, size_type col, difference_type nrows = -1, difference_type ncols = -1 ) const
332 {
333 difference_type nrows_ = rows( );
334 difference_type ncols_ = cols( );
335 if( nrows_ <= static_cast< difference_type >( row ) || nrows_ < static_cast< difference_type >( row + nrows ) )
336 {
337 return( false );
338 }
339 else if( ncols_ <= static_cast< difference_type >( col ) || ncols_ < static_cast< difference_type >( col + ncols ) )
340 {
341 return( false );
342 }
343 else if( nrows_ == nrows && ncols_ == ncols )
344 {
345 return( true );
346 }
347
348 if( nrows < 0 )
349 {
350 nrows = nrows_ - row;
351 }
352 if( ncols < 0 )
353 {
354 ncols = ncols_ - col;
355 }
356
357 if( out.resize( nrows, ncols ) )
358 {
359 const_pointer pi = paccess( row, col );
360 pointer po = out.paccess( 0, 0 );
361 for( difference_type c = 0 ; c < ncols ; c++ )
362 {
363 po = base::allocator_.copy_objects( pi, nrows, po );
364 pi += this->rows( );
365 }
366
367 return( true );
368 }
369 else
370 {
371 return( false );
372 }
373 }
374
385 bool trim( size_type row, size_type col, difference_type nrows = -1, difference_type ncols = -1 )
386 {
387 if( base::is_memory_shared( ) )
388 {
389 // 外部メモリを利用している場合
390 matrix o( *this );
391 return( o.trim( *this, row, col, nrows, ncols ) );
392 }
393 else
394 {
395 matrix o;
396
397 if( this->trim( o, row, col, nrows, ncols ) )
398 {
399 swap( o );
400
401 return( true );
402 }
403 else
404 {
405 return( false );
406 }
407 }
408 }
409
419 bool swap( matrix &m )
420 {
421 if( base::swap( m ) )
422 {
423 size_type _dmy_size1 = size1_;
424 size_type _dmy_size2 = size2_;
425 size1_ = m.size1_;
426 size2_ = m.size2_;
427 m.size1_ = _dmy_size1;
428 m.size2_ = _dmy_size2;
429 return( true );
430 }
431 else
432 {
433 return( false );
434 }
435 }
436
441 void clear( )
442 {
443 base::clear( );
444 size1_ = size2_ = 0;
445 }
446
447
448 size_type size1( ) const { return( size1_ ); }
449 size_type size2( ) const { return( size2_ ); }
450 size_type rows( ) const { return( size1_ ); }
451 size_type cols( ) const { return( size2_ ); }
452
453
454
455 /************************************************************************************************************
456 **
457 ** 行に対する順方向・逆方向の反復子
458 **
459 ************************************************************************************************************/
460
473 iterator row_begin( size_type r ){ return( iterator( paccess( r, 0 ), rows( ) ) ); }
474
476 const_iterator row_begin( size_type r ) const { return( const_iterator( paccess( r, 0 ), rows( ) ) ); }
477
479 iterator row_end( size_type r ){ return( iterator( paccess( r, cols( ) ), rows( ) ) ); }
480
482 const_iterator row_end( size_type r ) const { return( const_iterator( paccess( r, cols( ) ), rows( ) ) ); }
483
484
497 reverse_iterator row_rbegin( size_type r ){ return( reverse_iterator( row_end( r ) ) ); }
498
500 const_reverse_iterator row_rbegin( size_type r ) const { return( const_reverse_iterator( row_end( r ) ) ); }
501
503 reverse_iterator row_rend( size_type r ){ return( reverse_iterator( row_begin( r ) ) ); }
504
506 const_reverse_iterator row_rend( size_type r ) const { return( const_reverse_iterator( row_begin( r ) ) ); }
507
508
509 /************************************************************************************************************
510 **
511 ** 列に対する順方向・逆方向の反復子
512 **
513 ************************************************************************************************************/
514
527 iterator col_begin( size_type c ){ return( iterator( paccess( 0, c ), 1 ) ); }
528
530 const_iterator col_begin( size_type c ) const { return( const_iterator( paccess( 0, c ), 1 ) ); }
531
533 iterator col_end( size_type c ){ return( iterator( paccess( rows( ), c ), 1 ) ); }
534
536 const_iterator col_end( size_type c ) const { return( const_iterator( paccess( rows( ), c ), 1 ) ); }
537
538
551 reverse_iterator col_rbegin( size_type c ){ return( reverse_iterator( col_end( c ) ) ); }
552
554 const_reverse_iterator col_rbegin( size_type c ) const { return( const_reverse_iterator( col_end( c ) ) ); }
555
557 reverse_iterator col_rend( size_type c ){ return( reverse_iterator( col_begin( c ) ) ); }
558
560 const_reverse_iterator col_rend( size_type c ) const { return( const_reverse_iterator( col_begin( c ) ) ); }
561
562
563
564
565
566 public:
588 static const matrix diag( const value_type &s1, const value_type &s2, const value_type &s3 )
589 {
590 typedef __numeric__::value_compare< __numeric__::is_complex< value_type >::value > value_compare;
591 matrix d( 3, 3 );
592
593 if( value_compare::lt( s1, s2 ) )
594 {
595 // s1 < s2
596 if( value_compare::lt( s1, s3 ) )
597 {
598 if( value_compare::lt( s2, s3 ) )
599 {
600 d( 0, 0 ) = s3;
601 d( 1, 1 ) = s2;
602 d( 2, 2 ) = s1;
603 }
604 else
605 {
606 d( 0, 0 ) = s2;
607 d( 1, 1 ) = s3;
608 d( 2, 2 ) = s1;
609 }
610 }
611 else
612 {
613 d( 0, 0 ) = s2;
614 d( 1, 1 ) = s1;
615 d( 2, 2 ) = s3;
616 }
617 }
618 else
619 {
620 // s2 < s1
621 if( value_compare::lt( s1, s3 ) )
622 {
623 d( 0, 0 ) = s3;
624 d( 1, 1 ) = s1;
625 d( 2, 2 ) = s2;
626 }
627 else
628 {
629 if( value_compare::lt( s2, s3 ) )
630 {
631 d( 0, 0 ) = s1;
632 d( 1, 1 ) = s3;
633 d( 2, 2 ) = s2;
634 }
635 else
636 {
637 d( 0, 0 ) = s1;
638 d( 1, 1 ) = s2;
639 d( 2, 2 ) = s3;
640 }
641 }
642 }
643 return( d );
644 }
645
646
647
649 static const matrix _44(
650 const value_type &a00, const value_type &a01, const value_type &a02, const value_type &a03,
651 const value_type &a10, const value_type &a11, const value_type &a12, const value_type &a13,
652 const value_type &a20, const value_type &a21, const value_type &a22, const value_type &a23,
653 const value_type &a30, const value_type &a31, const value_type &a32, const value_type &a33
654 )
655 {
656 matrix o( 4, 4 );
657
658 o( 0, 0 ) = a00; o( 0, 1 ) = a01; o( 0, 2 ) = a02; o( 0, 3 ) = a03;
659 o( 1, 0 ) = a10; o( 1, 1 ) = a11; o( 1, 2 ) = a12; o( 1, 3 ) = a13;
660 o( 2, 0 ) = a20; o( 2, 1 ) = a21; o( 2, 2 ) = a22; o( 2, 3 ) = a23;
661 o( 3, 0 ) = a30; o( 3, 1 ) = a31; o( 3, 2 ) = a32; o( 3, 3 ) = a33;
662
663 return( o );
664 }
665
667 static const matrix _34(
668 const value_type &a00, const value_type &a01, const value_type &a02, const value_type &a03,
669 const value_type &a10, const value_type &a11, const value_type &a12, const value_type &a13,
670 const value_type &a20, const value_type &a21, const value_type &a22, const value_type &a23
671 )
672 {
673 matrix o( 3, 4 );
674
675 o( 0, 0 ) = a00; o( 0, 1 ) = a01; o( 0, 2 ) = a02; o( 0, 3 ) = a03;
676 o( 1, 0 ) = a10; o( 1, 1 ) = a11; o( 1, 2 ) = a12; o( 1, 3 ) = a13;
677 o( 2, 0 ) = a20; o( 2, 1 ) = a21; o( 2, 2 ) = a22; o( 2, 3 ) = a23;
678
679 return( o );
680 }
681
683 static const matrix _43(
684 const value_type &a00, const value_type &a01, const value_type &a02,
685 const value_type &a10, const value_type &a11, const value_type &a12,
686 const value_type &a20, const value_type &a21, const value_type &a22,
687 const value_type &a30, const value_type &a31, const value_type &a32
688 )
689 {
690 matrix o( 4, 3 );
691
692 o( 0, 0 ) = a00; o( 0, 1 ) = a01; o( 0, 2 ) = a02;
693 o( 1, 0 ) = a10; o( 1, 1 ) = a11; o( 1, 2 ) = a12;
694 o( 2, 0 ) = a20; o( 2, 1 ) = a21; o( 2, 2 ) = a22;
695 o( 3, 0 ) = a30; o( 3, 1 ) = a31; o( 3, 2 ) = a32;
696
697 return( o );
698 }
699
701 static const matrix _41( const value_type &a0, const value_type &a1, const value_type &a2, const value_type &a3 )
702 {
703 matrix o( 4, 1 );
704
705 o[ 0 ] = a0;
706 o[ 1 ] = a1;
707 o[ 2 ] = a2;
708 o[ 3 ] = a3;
709
710 return( o );
711 }
712
714 static const matrix _33(
715 const value_type &a00, const value_type &a01, const value_type &a02,
716 const value_type &a10, const value_type &a11, const value_type &a12,
717 const value_type &a20, const value_type &a21, const value_type &a22
718 )
719 {
720 matrix o( 3, 3 );
721
722 o( 0, 0 ) = a00; o( 0, 1 ) = a01; o( 0, 2 ) = a02;
723 o( 1, 0 ) = a10; o( 1, 1 ) = a11; o( 1, 2 ) = a12;
724 o( 2, 0 ) = a20; o( 2, 1 ) = a21; o( 2, 2 ) = a22;
725
726 return( o );
727 }
728
730 static const matrix _31( const value_type &a0, const value_type &a1, const value_type &a2 )
731 {
732 matrix o( 3, 1 );
733
734 o[ 0 ] = a0;
735 o[ 1 ] = a1;
736 o[ 2 ] = a2;
737
738 return( o );
739 }
740
742 static const matrix _22(
743 const value_type &a00, const value_type &a01,
744 const value_type &a10, const value_type &a11
745 )
746 {
747 matrix o( 2, 2 );
748
749 o( 0, 0 ) = a00; o( 0, 1 ) = a01;
750 o( 1, 0 ) = a10; o( 1, 1 ) = a11;
751
752 return( o );
753 }
754
756 static const matrix _21( const value_type &a0, const value_type &a1 )
757 {
758 matrix o( 2, 1 );
759
760 o[ 0 ] = a0;
761 o[ 1 ] = a1;
762
763 return( o );
764 }
765
767 static const matrix identity( size_type rows, size_type cols )
768 {
769 size_type size = rows < cols ? rows : cols;
770 matrix o( rows, cols );
771 for( size_type i = 0 ; i < size ; i++ )
772 {
773 o( i, i ) = 1;
774 }
775 return( o );
776 }
777
779 static const matrix zero( size_type rows, size_type cols )
780 {
781 return( matrix( rows, cols ) );
782 }
783
785 matrix operator -( ) const
786 {
787 const matrix &m = *this;
788 matrix o( size1_, size2_ );
789 for( size_type i = 0 ; i < o.size( ) ; i++ )
790 {
791 o[ i ] = -m[ i ];
792 }
793 return( o );
794 }
795
797 matrix t( ) const
798 {
799 const matrix &m = *this;
800 matrix o( size2_, size1_ );
801 for( size_type c = 0 ; c < o.cols( ) ; c++ )
802 {
803 for( size_type r = 0 ; r < o.rows( ) ; r++ )
804 {
805 o( r, c ) = m( c, r );
806 }
807 }
808 return( o );
809 }
810
812 matrix dagger( ) const
813 {
814 const matrix &m = *this;
815 matrix o( size2_, size1_ );
816 for( size_type c = 0 ; c < o.cols( ) ; c++ )
817 {
818 for( size_type r = 0 ; r < o.rows( ) ; r++ )
819 {
820 o( r, c ) = __numeric__::conjugate( m( c, r ) );
821 }
822 }
823 return( o );
824 }
825
826
827 // 行列に対する演算子
828 // += 行列
829 // += 定数
830 // -= 行列
831 // -= 定数
832 // *= 行列
833 // *= 定数
834 // /= 定数
835
844 template < class TT, class AAlocator >
845 const matrix& operator +=( const matrix< TT, AAlocator > &m2 )
846 {
847 matrix &m1 = *this;
848 #if _CHECK_ARRAY_OPERATION_ != 0
849 if( m1.size( ) != m2.size( ) )
850 {
851 // 足し算できません例外
852 ::std::cerr << "can't add arrays." << ::std::endl;
853 return( *this );
854 }
855 #endif
856 for( size_type i = 0 ; i < m1.size( ) ; i++ ) m1[i] += static_cast< T >( m2[i] );
857 return( m1 );
858 }
859
860
869 template < class TT, class AAlocator >
870 const matrix& operator -=( const matrix< TT, AAlocator > &m2 )
871 {
872 matrix &m1 = *this;
873 #ifdef _CHECK_ARRAY_OPERATION_
874 if( m1.size( ) != m2.size( ) )
875 {
876 // 引き算できません例外
877 ::std::cerr << "can't subtract matrixs." << ::std::endl;
878 return( m1 );
879 }
880 #endif
881 for( size_type i = 0 ; i < m1.size( ) ; i++ ) m1[i] -= static_cast< T >( m2[i] );
882 return( m1 );
883 }
884
885
894 template < class TT, class AAlocator >
895 const matrix& operator *=( const matrix< TT, AAlocator > &m2 )
896 {
897 matrix &m1 = *this;
898
899 #if defined( _CHECK_MATRIX_OPERATION_ ) && _CHECK_MATRIX_OPERATION_ != 0
900 if( m1.cols( ) != m2.rows( ) )
901 {
902 // 掛け算できません例外
903 ::std::cerr << "can't multiply matrices." << ::std::endl;
904 return( m1 );
905 }
906 #endif
907
908 typedef __numeric__::value_compare< __numeric__::is_complex< value_type >::value > value_compare;
909
910 matrix< T, Allocator > mat( m1.rows( ), m2.cols( ) );
911
912 #ifdef _OPENMP
913 int nCols = static_cast< int >( mat.cols( ) );
914
915 #pragma omp parallel for schedule( guided )
916 for( int c = 0 ; c < nCols ; c++ )
917 {
918 for( size_type t = 0 ; t < m1.cols( ) ; t++ )
919 {
920 value_type val = static_cast< value_type >( m2( t, c ) );
921 if( !value_compare::is_zero( val ) )
922 {
923 pointer pm0 = &mat( 0, c );
924 pointer pm1 = &m1( 0, t );
925 for( size_type r = 0 ; r < mat.rows( ) ; r++ )
926 {
927 pm0[ r ] += pm1[ r ] * val;
928 }
929 }
930 }
931 }
932 #else
933 size_type r, c, t;
934
935 for( c = 0 ; c < mat.cols( ) ; c++ )
936 {
937 for( t = 0 ; t < m1.cols( ) ; t++ )
938 {
939 value_type val = static_cast< value_type >( m2( t, c ) );
940 if( !value_compare::is_zero( val ) )
941 {
942 pointer pm0 = &mat( 0, c );
943 pointer pm1 = &m1( 0, t );
944 for( r = 0 ; r < mat.rows( ) ; r++ )
945 {
946 pm0[ r ] += pm1[ r ] * val;
947 }
948 }
949 }
950 }
951 #endif
952
953 m1.swap( mat );
954
955 return( m1 );
956 }
957
958
967 const matrix& operator +=( typename type_trait< T >::value_type val )
968 {
969 matrix &m = *this;
970 size_type i, size = m.rows( ) < m.cols( ) ? m.rows( ) : m.cols( );
971 for( i = 0 ; i < size ; i++ ) m( i, i ) += val;
972 return( m );
973 }
974
975
984 const matrix& operator -=( typename type_trait< T >::value_type val )
985 {
986 matrix &m = *this;
987 size_type i, size = m.rows( ) < m.cols( ) ? m.rows( ) : m.cols( );
988 for( i = 0 ; i < size ; i++ ) m( i, i ) -= val;
989 return( m );
990 }
991
992
1001 const matrix& operator *=( typename type_trait< T >::value_type val )
1002 {
1003 matrix &m = *this;
1004 for( size_type i = 0 ; i < m.size( ) ; i++ ) m[i] *= val;
1005 return( m );
1006 }
1007
1008
1017 const matrix& operator /=( typename type_trait< T >::value_type val )
1018 {
1019 matrix &m = *this;
1020 #if defined( _CHECK_MATRIX_OPERATION_ ) && _CHECK_MATRIX_OPERATION_ != 0
1021 if( val == value_type( 0 ) )
1022 {
1023 // ゼロ除算発生
1024 ::std::cerr << "zero division occured." << ::std::endl;
1025 return( m );
1026 }
1027 #endif
1028 for( size_type i = 0 ; i < m.size( ) ; i++ ) m[i] /= val;
1029 return( m );
1030 }
1031
1032
1041 bool is_equal( const matrix &a, const double delta )
1042 {
1043 typedef __numeric__::value_compare< __numeric__::is_complex< value_type >::value > value_compare;
1044
1045 if( rows( ) != a.rows( ) || cols( ) != a.cols( ) )
1046 {
1047 return( false );
1048 }
1049
1050 for( size_type i = 0 ; i < base::size( ) ; i++ )
1051 {
1052 if( !value_compare::eq( operator []( i ), a[ i ], delta ) )
1053 {
1054 return( false );
1055 }
1056 }
1057
1058 return( true );
1059 }
1060
1068 bool operator ==( const matrix &a ) const
1069 {
1070 if( rows( ) != a.rows( ) || cols( ) != a.cols( ) )
1071 {
1072 return( false );
1073 }
1074
1075 for( size_type i = 0 ; i < base::size( ) ; i++ )
1076 {
1077 if( operator []( i ) != a[ i ] )
1078 {
1079 return( false );
1080 }
1081 }
1082
1083 return( true );
1084 }
1085
1093 bool operator !=( const matrix &a ) const { return( !( *this == a ) ); }
1094
1095
1096 public:
1107 template < class TT, class AAlocator >
1108 const matrix& operator =( const matrix< TT, AAlocator > &o )
1109 {
1110 base::operator =( o );
1111
1112 if( base::empty( ) )
1113 {
1114 size1_ = size2_ = 0;
1115 }
1116 else
1117 {
1118 size1_ = o.size1( );
1119 size2_ = o.size2( );
1120 }
1121
1122 return( *this );
1123 }
1124
1134 const matrix< T, Allocator >& operator =( const matrix< T, Allocator > &o )
1135 {
1136 if( this == &o ) return( *this );
1137
1138 base::operator =( o );
1139
1140 if( base::empty( ) )
1141 {
1142 size1_ = size2_ = 0;
1143 }
1144 else
1145 {
1146 size1_ = o.size1( );
1147 size2_ = o.size2( );
1148 }
1149
1150 return( *this );
1151 }
1152
1153
1154 // 要素へのアクセス
1155 private:
1163 pointer paccess( size_type r, size_type c )
1164 {
1165 return( base::data_ + r + c * size1_ );
1166 }
1167
1175 const_pointer paccess( size_type r, size_type c ) const
1176 {
1177 return( base::data_ + r + c * size1_ );
1178 }
1179
1180
1181 public:
1191 reference at( size_type r, size_type c )
1192 {
1193 _CHECK_ACCESS_VIOLATION2U_( r, c )
1194 return( base::data_[ r + c * size1_ ] );
1195 }
1196
1206 const_reference at( size_type r, size_type c ) const
1207 {
1208 _CHECK_ACCESS_VIOLATION2U_( r, c )
1209 return( base::data_[ r + c * size1_ ] );
1210 }
1211
1221 reference operator ()( size_type r, size_type c )
1222 {
1223 _CHECK_ACCESS_VIOLATION2U_( r, c )
1224 return( base::data_[ r + c * size1_ ] );
1225 }
1226
1236 const_reference operator ()( size_type r, size_type c ) const
1237 {
1238 _CHECK_ACCESS_VIOLATION2U_( r, c )
1239 return( base::data_[ r + c * size1_ ] );
1240 }
1241
1242
1243 public:
1245 matrix( ) : base( ), size2_( 0 ), size1_( 0 ) {}
1246
1248 explicit matrix( const Allocator &a ) : base( a ), size2_( 0 ), size1_( 0 ) {}
1249
1250
1252 matrix( size_type rnum, size_type cnum ) : base( rnum * cnum ), size2_( cnum ), size1_( rnum )
1253 {
1254 if( base::empty( ) ) size1_ = size2_ = 0;
1255 }
1256
1258 matrix( size_type rnum, size_type cnum, const Allocator &a ) : base( rnum * cnum, a ), size2_( cnum ), size1_( rnum )
1259 {
1260 if( base::empty( ) ) size1_ = size2_ = 0;
1261 }
1262
1263
1265 matrix( size_type rnum, size_type cnum, const T &val ) : base( rnum * cnum, val ), size2_( cnum ), size1_( rnum )
1266 {
1267 if( base::empty( ) ) size1_ = size2_ = 0;
1268 }
1269
1271 matrix( size_type rnum, size_type cnum, const T &val, const Allocator &a ) : base( rnum * cnum, val, a ), size2_( cnum ), size1_( rnum )
1272 {
1273 if( base::empty( ) ) size1_ = size2_ = 0;
1274 }
1275
1276
1277
1278
1280 matrix( size_type rnum, size_type cnum, pointer ptr, size_type mem_available ) : base( rnum * cnum, ptr, mem_available ), size2_( cnum ), size1_( rnum )
1281 {
1282 if( base::empty( ) ) size1_ = size2_ = 0;
1283 }
1284
1286 matrix( size_type rnum, size_type cnum, const T &val, pointer ptr, size_type mem_available ) : base( rnum * cnum, val, ptr, mem_available ), size2_( cnum ), size1_( rnum )
1287 {
1288 if( base::empty( ) ) size1_ = size2_ = 0;
1289 }
1290
1291
1296 template < class TT, class AAlocator >
1297 matrix( const matrix< TT, AAlocator > &o ) : base( o ), size2_( o.size2( ) ), size1_( o.size1( ) )
1298 {
1299 if( base::empty( ) ) size1_ = size2_ = 0;
1300 }
1301
1302
1304 matrix( const matrix< T, Allocator > &o ) : base( o ), size2_( o.size2_ ), size1_( o.size1_ )
1305 {
1306 if( base::empty( ) ) size1_ = size2_ = 0;
1307 }
1308 };
1309
1310
1311
1325 template < class T, class Allocator >
1326 inline ::std::ostream &operator <<( ::std::ostream &out, const matrix< T, Allocator > &m )
1327 {
1328 typename matrix< T, Allocator >::size_type r, c;
1329 for( r = 0 ; r < m.rows( ) ; r++ )
1330 {
1331 if( r != 0 )
1332 {
1333 out << ::std::endl;
1334 }
1335 for( c = 0 ; c < m.cols( ) ; c++ )
1336 {
1337 out << m( r, c );
1338 if( c != m.cols( ) - 1 ) out << ", ";
1339 }
1340 }
1341
1342 return( out );
1343 }
1344
1345
1346 /************************************************************************************************************
1347 **
1348 ** 行列に対する演算子
1349 ** 行列 + 行列
1350 ** 行列 + 定数
1351 ** 定数 + 行列
1352 **
1353 ** 行列 - 行列
1354 ** 行列 - 定数
1355 ** 定数 - 行列
1356 **
1357 ** 行列 * 行列
1358 ** 行列 * 定数
1359 ** 定数 * 行列
1360 **
1361 ** 行列 / 定数
1362 **
1363 ************************************************************************************************************/
1364
1365
1373 template < class T, class Allocator >
1374 inline matrix< T, Allocator > operator +( const matrix< T, Allocator > &m1, const matrix< T, Allocator > &m2 )
1375 {
1376 return( matrix< T, Allocator >( m1 ) += m2 );
1377 }
1378
1379
1387 template < class T, class Allocator >
1388 inline matrix< T, Allocator > operator -( const matrix< T, Allocator > &m1, const matrix< T, Allocator > &m2 )
1389 {
1390 return( matrix< T, Allocator >( m1 ) -= m2 );
1391 }
1392
1393
1401 template < class T, class Allocator >
1402 inline matrix< T, Allocator > operator *( const matrix< T, Allocator > &m1, const matrix< T, Allocator > &m2 )
1403 {
1404 return( matrix< T, Allocator >( m1 ) *= m2 );
1405 }
1406
1407
1417 template < class T, class Allocator >
1418 inline matrix< T, Allocator > operator +( const matrix< T, Allocator > &m, typename type_trait< T >::value_type val )
1419 {
1420 return( matrix< T, Allocator >( m ) += val );
1421 }
1422
1432 template < class T, class Allocator >
1433 inline matrix< T, Allocator > operator +( typename matrix< T, Allocator >::value_type val, const matrix< T, Allocator > &m )
1434 {
1435 return( matrix< T, Allocator >( m ) += val );
1436 }
1437
1438
1439
1449 template < class T, class Allocator >
1450 inline matrix< T, Allocator > operator -( const matrix< T, Allocator > &m, typename type_trait< T >::value_type val )
1451 {
1452 return( matrix< T, Allocator >( m ) -= val );
1453 }
1454
1464 template < class T, class Allocator >
1465 inline matrix< T, Allocator > operator -( typename matrix< T, Allocator >::value_type val, const matrix< T, Allocator > &m )
1466 {
1467 return( matrix< T, Allocator >( m ) -= val );
1468 }
1469
1470
1480 template < class T, class Allocator >
1481 inline matrix< T, Allocator > operator *( const matrix< T, Allocator > &m, typename type_trait< T >::value_type val )
1482 {
1483 return( matrix< T, Allocator >( m ) *= val );
1484 }
1485
1495 template < class T, class Allocator >
1496 inline matrix< T, Allocator > operator *( typename matrix< T, Allocator >::value_type val, const matrix< T, Allocator > &m )
1497 {
1498 return( matrix< T, Allocator >( m ) *= val );
1499 }
1500
1501
1502 // @brief 定数との割り算
1511 template < class T, class Allocator >
1512 inline matrix< T, Allocator > operator /( const matrix< T, Allocator > &m, typename type_trait< T >::value_type val )
1513 {
1514 return( matrix< T, Allocator >( m ) /= val );
1515 }
1516
1517
1518 // mist名前空間の終わり
1519 _MIST_END
1520
1521 #endif // __INCLUDE_MIST_MATRIX__
1522