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 
Generated on Wed Nov 12 2014 19:44:18 for MIST by doxygen 1.8.1.2